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Abstract 

A  total-field  scattered-field  (TFSF)  boundary  can  be  used  to  introduce  incident  plane  waves 
into  finite-difference  time-domain  (FDTD)  simulations.  For  fields  which  are  traveling  obliquely 
to  the  grid  axes,  there  is  no  simple  way  to  account  fully  for  the  effects  of  the  inherent  numeri¬ 
cal  artifacts  associated  with  plane-wave  propagation  in  the  FDTD  grid.  Failure  to  account  for 
these  artifacts  causes  erroneous  fields  to  leak  across  the  TFSF  boundary.  The  work  performed 
under  this  grant  developed  a  TFSF  boundary  for  problems  involving  a  planar  interface  that  is 
theoretically  exact.  Such  a  boundary  is  essential  when  modeling  the  scattering  from  objects 
in  the  vicinity  of  ocean  sediment.  Although  the  TFSF  boundary  assumes  a  planar  interface, 
the  actual  interface  used  in  a  simulation  can  take  any  shape  and  the  space  within  the  total-field 
region  is  arbitrary.  A  suite  of  C  programs  was  written  to  implement  the  boundary. 


1  Introduction 


The  total-field  scattered-field  (TFSF)  formulation  is  a  method  for  introducing  energy  into  a  finite- 
difference  time-domain  (FDTD)  simulation.  It  defines  a  boundary,  identified  as  the  TFSF  bound¬ 
ary,  which  divides  the  computational  domain  into  two  regions:  a  total-field  (TF)  region  which 
contains  both  the  incident  field  and  any  scattered  field,  and  a  scattered-field  (SF)  region  which 
contains  only  scattered  fields.  Scatterers  are  confined  to  exist  within  the  TF  region.  Throughout 
the  grid  the  fields  must  have  self-consistent  update  equations  meaning  that  nodes  in  the  TF  region 
must  depend  on  the  total  field  at  neighboring  nodes  while  nodes  in  the  SF  region  must  depend  on 
the  scattered  field  at  neighboring  nodes.  However,  nodes  which  are  tangential  to  the  TFSF  bound¬ 
ary  will  have  at  least  one  neighboring  node  in  the  region  different  from  their  own.  Rectification  of 
this  inconsistency  is  what  drives  the  TFSF  method. 

Given  knowledge  of  the  incident  field  which  should  exist  at  nodes  tangential  to  the  TFSF 
boundary,  one  does  the  following.  For  the  update  of  a  node  which  is  in  the  TF  region  and  depends 
on  a  neighboring  node  in  the  SF  region,  the  incident  field  is  added  to  that  neighboring  node. 
Conversely,  for  the  update  of  a  node  which  is  in  the  SF  region  and  depends  on  a  neighboring 
node  in  TF  region,  the  incident  field  is  subtracted  from  that  neighboring  node.  In  this  way  the 
TFSF  boundary  acts  as  a  Huygens  surface  and  was  originally  described  as  such  by  Merewether  et 
al  [1]. 
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In  order  to  implement  the  TFSF  method,  one  must  know  the  incident  field  at  every  time-step 
at  all  the  tangential  nodes  adjacent  to  the  boundary.  Historically  the  TFSF  method  has  been  used 
to  introduce  plane  waves  into  the  grid,  but  in  theory  any  incident  field  could  be  realized.  In  this 
work  we  restricted  ourselves  to  plane  waves.  In  the  continuous  world,  the  analytic  description  of 
pulsed  plane-wave  propagation  is  relatively  trivial.  Unfortunately,  because  of  the  inherent  differ¬ 
ence  between  the  way  in  which  waves  propagate  in  the  FDTD  grid  and  the  way  in  which  they 
propagate  in  the  continuous  world,  one  should  not  simply  use  the  continuous-world  expression  for 
the  incident  plane  wave.  If  one  were  to  do  that,  the  mismatch  between  the  analytic  description 
of  the  incident  field  and  how  the  incident  field  actually  propagates  within  the  grid  causes  fields  to 
leak  across  the  TFSF  boundary  (in  the  absence  of  a  scatterer,  no  fields  should  be  present  in  the  SF 
region).  Fortunately  there  is  a  relatively  simple  fix  to  this  problem. 

An  auxiliary  one-dimensional  (ID)  FDTD  simulation  can  be  implemented  to  model  the  prop¬ 
agation  of  the  incident  plane  wave.  If  the  propagation  direction  of  the  incident  field  in  a  higher¬ 
dimensional  grid  is  aligned  with  one  of  the  grid  axes,  the  auxiliary  ID  grid  can  be  used  to  describe 
exactly  the  incident  field  at  all  nodes  adjacent  to  the  TFSF  boundary.  There  is  no  leakage.  A 
detailed  discussion  of  the  implementation  can  be  found  in  [2]. 

When  the  incident  field  propagates  obliquely,  a  one-dimensional  auxiliary  grid  can  still  be  used 
to  describe,  at  least  approximately,  the  incident  field.  However,  using  these  fields  in  the  higher¬ 
dimensional  simulation  has  inherent  errors.  First,  interpolation  must  be  used  to  find  the  fields  at 
points  on  the  ID  grid  corresponding  to  the  projected  locations  of  nodes  in  the  higher-dimensional 
grid.  Implementation  details  for  oblique  propagation  can  also  be  found  in  [2].  Second,  the  disper¬ 
sion  which  the  fields  experience  in  the  ID  grid  does  not  correspond  to  the  dispersion  experienced  in 
the  higher-dimensional  grid.  To  help  rectify  this,  Guiffaut  and  Mahdjoubi  [3]  proposed  a  technique 
which  modified  the  Courant  number  in  the  auxiliary  grid  so  that  the  dispersion  nearly  matched  that 
of  the  higher-dimensional  grid  (see  also  Sec.  5.9.1  of  [2]).  Nevertheless,  there  are  still  differences 
between  propagation  in  the  two  grids  and  the  need  for  interpolation,  which  inherently  introduces 
errors,  still  exists.  Third,  the  orientation  of  vector  fields  (i.e.,  the  velocity  field  in  acoustic  simula¬ 
tions  and  either  the  electric  of  magnetic  field  in  electromagnetic  simulations)  is  not  the  same  in  the 
FDTD  grid  as  it  is  in  the  continuous  world.  This  fact  is  discussed  in  [4,5]  and  had  been  considered 
previously  in  [6, 7].  The  orientation  is  dependent  on  frequency  and  hence  a  simple  scalar  cannot 
be  used  to  project  field  components  from  a  ID  grid  to  a  higher-dimensional  grid. 

Instead  of  employing  an  auxiliary  grid,  two  recent  papers  proposed  a  technique  where  the 
incident  field  is  obtained  analytically  by  way  of  the  FDTD  dispersion  relation  [4,5].  To  distinguish 
this  approach  from  that  which  relies  upon  auxiliary  grids,  we  label  this  approach  the  analytic  field 
propagation  (AFP)  TFSF  technique.  Moss  et  al.  considered  the  situation  of  an  electromagnetic 
field  incident  on  layered  uniaxial  anisotropic  media  and  hence  had  to  account  for  the  reflection  and 
transmission  coefficients  of  the  layers  [4].  Schneider  restricted  consideration  to  propagation  in  a 
homogeneous  space  [5]. 

Because  of  the  complexity  of  the  media  being  considered  by  Moss  et  al.,  the  equations  given 
in  their  work  masked  the  simplicity  which  pertains  for  problems  involving  isotropic  media  and 
problems  involving  a  plane  satisfying  a  Dirichlet  boundary  condition.  A  Dirichlet  boundary  in 
an  electromagnetic  simulation  is  equivalent  to  either  a  pressure-release  or  rigid  boundary  in  an 
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acoustic  simulation  depending  on  the  polarization  assumed.  We  identify  such  a  boundary  a  perfect 
electrically  conducting  (PEC)  boundary.  In  this  report  we  examined  these  problems  in  some  detail 
and  obtained  equations  which  are  significantly  simpler  than  those  presented  for  anisotropic  media. 
As  described  both  by  Moss  et  al.  [4]  and  Schneider  [5],  the  AFP  TFSF  method  requires  Fourier 
transforms  to  obtain  the  incident  field.  However,  neither  of  those  papers  provided  details  of  the 
type  of  transform  actually  need.  In  fact,  an  exact  implementation  requires  not  a  discrete  Fourier 
transform  but  rather  a  discrete-time  Fourier  transform  (which  involves  a  continuous  integral).  As 
will  be  discussed,  this  transform  can  be  approximated  by  a  discrete  transform,  but  one  should  be 
aware  of  the  inherent  approximation. 

This  report  begins  by  discussing  the  relationship  between  acoustic  and  electromagnetic  FDTD 
simulations.  The  remainder  of  the  work  is  presented  in  terms  of  electromagnetics.  An  overview 
of  the  AFP  TFSF  technique  is  then  presented  that  discusses  the  issues  concerning  the  Fourier 
transform.  We  then  consider  problems  involving  a  PEC  plane  or  a  boundary  between  two  half 
spaces.  Both  TE  and  TM  polarization  are  considered.  As  will  be  shown,  these  two-dimensional 
simulations  are  analogs  of  the  acoustic  problem  where  the  single  scalar  field  (i.e.,  the  magnetic 
field  in  the  case  of  TE  polarization  and  electric  field  in  the  case  of  TM  polarization)  can  be  equated 
with  pressure  and  the  vector  field  with  velocity.  Initially  we  restrict  ourselves  to  lossless  material 
but  then  consider  lossy  material  and  incidence  angles  beyond  the  critical  angle.  Finally,  there  is 
a  discussion  of  the  use  of  the  the  AFP  TFSF  technique  in  three  dimensional  simulations.  Part  of 
the  work  described  here  has  been  accepted  for  publication  in  the  IEEE  Transactions  on  Antennas 
and  Propagation.  Another  paper,  also  based  on  the  work  presented  here,  has  been  submitted  to  the 
same  journal  and  is  currently  in  review. 


2  Relationship  between  Acoustic  and  Electromagnetic  FDTD 
Grids 


The  FDTD  method  employs  finite-differences  to  approximate  Ampere’s  and  Faraday’s  laws.  Am¬ 
pere’s  and  Faraday’s  laws  are  first-order  differential  equations  which  couple  the  electric  and  mag¬ 
netics  fields.  As  we  have  seen,  with  a  judicious  discretization  of  space  and  time,  the  resulting 
equations  can  be  solved  for  “future”  fields  in  terms  of  known  past  fields. 

Other  physical  phenomena  are  also  described  by  coupled  first-order  differential  equations 
where  the  temporal  derivative  of  one  field  is  related  to  the  spatial  derivative  of  another  field.  Both 
acoustics  and  elastic  wave  propagation  are  such  phenomena.  Here  we  will  consider  only  acous¬ 
tic  propagation.  Specifically  we  will  consider  small-signal  acoustics  which  can  be  described  in 
terms  of  the  scalar  pressure  field  P(x ,  y.  z,  t)  and  the  vector  velocity  v(x,  y,  z,  t ).  The  material 
parameters  are  the  speed  of  sound  ca  and  the  density  p  (both  of  which  can  vary  as  a  function  of 
position). 

The  governing  acoustic  equations  in  three  dimensions  are 


8P 
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dt 
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or,  expanded  in  terms  of  the  components, 
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Equation  (2)  is  essentially  a  variation  of  Newton’s  second  law,  F  =  ma,  where  instead  of  accel¬ 
eration  a  there  is  the  derivative  of  the  velocity,  instead  of  mass  m  there  is  the  mass  density,  and 
instead  of  force  F  there  is  the  derivative  of  pressure  (which  is  a  force  per  area  and  the  negative 
sign  accounts  for  the  fact  that  if  pressure  is  building  in  a  particular  direction  that  tends  to  cause  ac¬ 
celeration  in  the  opposite  direction).  Equation  (1)  comes  from  an  equation  of  state  for  the  material 
(with  various  approximations  assumed  along  the  way). 

Taking  the  divergence  of  (2)  and  interchanging  the  order  of  temporal  and  spatial  differentiation 
yields 

|-V.v=--V2P.  (7) 

at  p 


Taking  the  temporal  derivative  of  (1)  and  using  (7)  yields 

d2P  n  9  ,, 

li 5-  =  -^SV-v  =  c«V-P'  <8> 

Rearranging  this  yields  the  wave  equation 

V2P  -  — —  =  o.  (9) 

cl  dt 2  w 

Thus  the  usual  techniques  and  solutions  one  is  familiar  with  from  electromagnetics  carry  over  to 
acoustics.  For  example,  a  harmonic  plane  wave  given  by 

P(x,  y,  z,  t)  —  Poe-jk'reJW*  (10) 

is  a  valid  solution  to  the  governing  equations  where  P0  is  a  constant  and  the  wave  vector  k  can  be 
written 

k  =  kxax  +  kySLy  +  kz  az  =  kak  =  (u;/ca)  &k.  (11) 

Substituting  (10)  into  (4)  and  assuming  exp(jut)  temporal  dependence  yields 


Rearranging  terms  gives 


jujvx  =  ~{—jkx)P. 
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Following  the  same  steps  for  the  y  and  z  components  produces 


Vy  = 


=  ^ P, 


r, 

V,  =  - P. 


Thus  the  harmonic  velocity  is  given  by 

1  k 

v  =  vxax  +  vvay  +  vzaz  =  — (kx  ax  +  kyay  +  kzaz)P  =  — Pak.  (16) 

PU)  pu! 

Since  the  wave  number,  i.e.,  the  magnitude  of  the  wave  vector,  is  given  by  k  —  Lo/ca,  the  ratio  of 
the  magnitude  of  pressure  to  the  velocity  is  given  by 

~  =PCa-  (17) 

The  term  on  the  right-hand  side  is  known  as  the  characteristic  impedance  of  the  medium  which  is 
often  written  as  Z. 


2.1  Governing  FDTD  Equations 

To  obtain  an  FDTD  algorithm  for  acoustic  propagation,  the  pressure  and  components  of  velocity 
are  discretized  in  both  time  and  space.  In  electromagnetics  there  are  two  vector  fields  and  hence 
six  field-components  which  have  to  be  arranged  in  space-time.  In  acoustics  there  is  one  scalar  field 
and  one  vector  field.  Thus  there  are  only  four  field-components. 

To  implement  a  3D  acoustic  FDTD  algorithm,  a  suitable  arrangement  of  nodes  is  as  shown 
in  Fig.  1.  A  pressure  node  is  surrounded  by  velocity  components  such  that  the  components  are 
oriented  along  the  line  joining  the  component  and  the  pressure  node.  This  should  be  contrasted 
to  the  arrangement  of  nodes  in  electromagnetic  grids  where  the  components  of  the  magnetic  field 
swirled  around  the  components  of  the  electric  field,  and  vice  versa.  In  electromagnetics  one  is 
modeling  coupled  curl  equations  where  the  partial  derivatives  are  related  to  behavior  orthogonal  to 
the  direction  of  the  derivative.  In  acoustics,  where  the  governing  equations  involve  the  divergence 
and  gradient,  the  partial  derivatives  are  associated  with  behavior  in  the  direction  of  the  derivative. 

The  arrangement  of  nodes  in  a  2D  grid  is  illustrated  in  Fig.  2.  (This  should  can  be  compared 
to  the  2D  electromagnetic  grids  which  will  be  discussed  later,  i.e.,  Fig.  6  for  the  TEZ  case  and  Fig. 
8  for  the  TMZ  case.)  Because  pressure  is  inherently  an  acoustic  field,  there  are  not  two  different 
polarization  associated  with  2D  acoustic  simulations — nor  is  there  a  notion  of  polarization  in  three 
dimensions. 

In  addition  to  the  spatial  offsets,  the  pressure  nodes  are  assumed  to  be  offset  a  half  temporal 
step  from  the  velocity  nodes  (but  all  the  velocity  components  exist  at  the  same  time-step).  The 
following  notation  will  be  used  with  an  implicit  understanding  of  spatial  offsets 

p(x,y,z,t)  =  P(mAx,nAy,pAz,qAt)  =  Pq[m,n,p],  (18) 

vx(x,y,  z,t)  =  vx([m  +  l/2]Ax,nAy,pAz,[q+l/2]At)  =  vq+l/2[m,n,p\,  (19) 

Vy(x,y,z,t)  =  Vy  (mAx,  [n— t—  l/2]Ay,pAz,  [q  +  1/2} At)  =  vq+1/2[m,n,p],  (20) 

vz(x,y,z,t)  =  vz(mAx,nAy,\p  +  l/2]Az,[q  +  l/2]At)  =  vqz+1/2[m,n,p\.  (21) 
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Figure  1:  An  acoustic  unit  cell  in  three  dimensions  showing  the  arrangement  of  velocity  nodes 
relative  to  the  pressure  node  with  the  same  spatial  indices. 


We  will  assume  the  spatial  step  sizes  are  the  same,  i.e.,  Ax  =  Ay  —  Az  =  5. 

Replacing  the  derivatives  in  (3)  with  finite  differences  and  using  the  discretization  of  (18)— (21) 
yields  the  following  update  equation: 

Pg[m,n,p }  =  Pq^[m,n,p]- pcl^  (vgx-1/2[m,n,p]-vqx-1/2[m-l,n,p]  + 

vq-1/2[m,n,p\  -  vq-1/2[m,n  -  1  ,p]  + 
vq-1/2[m,n,p]  -  v9z-1/2[m,n,p  -  1]) .  (22) 


The  sound  speed  and  the  density  can  be  functions  of  space.  Let  us  assume  that  the  density  and 
sound  speed  are  specified  at  the  grid  points  corresponding  to  the  location  of  pressure  nodes.  Addi¬ 
tionally,  assume  that  the  sound  speed  can  be  defined  in  terms  of  a  background  sound  speed  c0  and 
a  relative  sound  speed  cy: 

ca  —  rycQ.  (23) 

The  background  sound  speed  corresponds  to  the  fastest  speed  of  propagation  at  any  location  in  the 
grid  so  that  <y  <  1.  The  coefficient  of  the  spatial  finite-difference  in  (22)  can  now  be  written 


2  A t  2  Co  At  2  a 

pca—  =  pc;c0-j-  =  pcrcoSc 


(24) 


where,  similar  to  electromagnetics,  the  Courant  number  is  Sc  =  c0At/8.  The  explicit  spatial 
dependence  of  the  density  and  sound  speed  can  be  emphasized  by  writing  the  coefficient  as 


p[m,  n,  p]c2  [m,  n,  p]c0Sc 


(25) 
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Figure  2:  The  arrangement  of  nodes  in  a  2D  acoustic  simulation.  In  a  computer  FDTD  implemen¬ 
tation  the  nodes  shown  within  the  dashed  enclosures  will  have  the  same  spatial  indices.  This  is 
illustrated  by  the  two  depictions  of  a  unit  cell  at  the  bottom  of  the  figure.  The  one  on  the  left  shows 
the  nodes  with  the  spatial  offsets  given  explicitly.  The  one  on  the  right  shows  the  corresponding 
node  designations  which  would  be  used  in  a  computer  program.  (Here  Pr  is  used  for  the  pressure 
array.) 
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where  p[m,n,p\  is  the  density  which  exists  at  the  same  point  as  the  pressure  node  P[m,n,p) 
and  Cr[m,  n,p ]  is  the  relative  sound  speed  which  exists  at  this  same  point.  Note  that  the  Courant 
number  Sc  and  the  background  sound  speed  c0  are  independent  of  position.  Furthermore,  the  entire 
coefficient  is  independent  of  time. 

The  update  equation  for  the  x  component  of  velocity  is  obtained  from  the  discretized  version 
of  (4)  which  yields 

vq+1/2[m,n,p]  =  vqx-l/2[m,n,p]  -  (Pq[m  +  l,n,p\  -  Pq[m,n,p ])  (26) 

The  coefficient  of  this  equation  does  not  contain  the  Courant  number  but  that  can  be  obtained  by 
multiplying  and  dividing  by  the  background  sound  speed 


1  At  _  1  c0  At  _  J_ 
p  S  pco  5  pco 


(27) 


We  wish  to  define  the  density  only  at  the  pressure  nodes.  Since  the  x -component  of  the  velocity 
is  offset  from  the  pressure  a  half  spatial  step  in  the  x  direction,  what  is  the  appropriate  velocity 
to  use?  The  answer,  much  as  it  is  in  the  case  of  an  interface  between  two  different  materials  in 
electromagnetics,  is  the  average  of  the  densities  to  either  side  of  the  pressure  node  (where  the 
notion  of  “either  side”  is  dictated  by  the  orientation  of  the  velocity  node).  Therefore  the  coefficient 
can  be  written 


_ 2 Sc _ 

(p[m  +  l,n,p]  + p[m,n,p])c0' 


(28) 


The  update  equations  for  the  velocity  components  can  now  be  written  as 


vq+1/2[m,n,p] 
vq+1/2[m,  n,p] 
vq+1/2[m,n,p] 


vq  1/2  [m,  n,  p}- 

2  s 

7-7 - x - A - - - — —  (Pq[m  +  1,  n,p]  —  Pq[m,  n. pi) , 

(p[m,n,p}+p[m+l,n,p})coK  1  ’  1 

vq~l/2[m,n,p]  - 

2  s 

(p[m,n,p\+  p[m,n+l,p])c0K  J; 

vq~1/2[m,n,p]  - 

2  s 

7-7 - x - x~ - -rr —  (Pq[m,n,p  +  1]  -  Pq[m,n,p]) . 

(p[m,n,p]  +p[m,n,p+  l])c0  v 


(29) 

(30) 

(31) 


2.2  Two-Dimensional  Implementation 

Let  us  consider  a  2D  simulation  in  which  the  fields  vary  in  the  x  and  y  directions.  The  grid  would 
be  as  shown  in  Fig.  2  and  it  is  assumed  that  Ax  =  A y  —  S.  Assume  the  arrays  pr,  vx,  and  vy 
hold  the  pressure,  x  component  of  the  velocity,  and  the  y  component  of  the  velocity,  respectively. 
Assume  the  macros  Pr,  Vx,  and  Vy  have  been  created  to  facilitate  accessing  these  arrays.  The 
update  equations  can  be  written 
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Vx(m,n)  =  Vx(m,n)  -  Cvxp (m, n) * {Pr (m+1 , n)  -  Pr(m,n)); 
Vy(m,n)  =  Vy(m,n)  -  Cvyp (m, n) * (Pr (m, n+1 )  -  Pr(m,n) ) ; 
Pr(m,n)  =  Pr(m,n)  -  Cprv  (m,  n)  *  {  (Vx  (m,  n)  -  Vx(m-l,n)) 

+  {Vy (m,  n)  -  Vy(m,n-1))) 


where  the  coefficient  arrays  are  given  by 

C vxpix^t  n.)  —  Sc 


±Sc  2Sc 

PC 0  '  {m+l/2)S,nS  (p[™  +!,«]+  PK  «])cq  ’ 


1  2  S 

Cvyp(m,  n)  =  Sc  =  r— .  ;  :  f  .. -i\  >  (33) 

PCo  m8, (n+l/2)5  (p[m,  Vl\  +  p[m,  U  +  l])c0 

Cprv(m,n )  =  p[m,n]Cr[m)n]co5'c.  (34) 

These  update  equations  are  little  different  from  those  for  the  TMZ  case.  The  TMZ  update 
equations  are 

Hy(m,n)  =  Hy(m,n)  +  Chye (m, n) * (Ez (m+1 , n)  -  Ez(m,n) ) ; 

Hx(m,n)  =  Hx(m,n)  -  Chxe (m,n) * (Ez (m,n+l)  -  Ez(m,n)); 

Ez(m,n)  =  Ez(m,n)  +  Cezh (m, n) * ( (Hy (m, n)  -  Hy(m-l,n)) 

-(Hx(m,n)  -  Hx (m, n-1) ) ) ; 

There  is  a  one-to-one  mapping  between  these  sets  of  equations.  One  can  equate  values  as  follows 


vx 

~Hy> 

(35) 

Vy 

<=> 

Hx, 

(36) 

P 

Ez, 

(37) 

Cvxp 

Chye i 

(38) 

Cyyp 

Chxe, 

(39) 

Cprv 

Cezh ■ 

(40) 

Thus,  converting  2D  programs  which  were  written  to  model  electromagnetic  field  propagation  to 
ones  which  can  model  acoustic  propagation  is  surprisingly  straightforward.  Essentially,  all  one  has 
to  do  is  change  some  labels  and  a  few  signs. 

For  TEZ  simulations,  the  updated  equations  for  a  lossless  medium  are 

Hz(m,n)  =  Hz(m,n)  + 

Chze(m,n) * ( (Ex(m,n+1) -Ex(m,n) ) - (Ey(m+l,n) -Ey(m,n) ) ) ; 
Ex(m,n)  =  Ex(m,n)  +  Cexh (m, n) * (Hz (m, n) -Hz (m, n-1) ) ; 

Ey(m,n)  =  Ey(m,n)  -  Ceyh (m, n) * (Hz (m; n) -Hz (m-1 , n) ) ; 
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In  this  case  the  conversion  from  the  electromagnetic  equations  to  the  acoustic  equations  can  be 
accomplished  with  the  following  mapping 


Ey, 

(41) 

Vy 

~EX, 

(42) 

P 

<£> 

Hz, 

(43) 

Cvxp 

Ceyhi 

(44) 

F vyp 

E  exhi 

(45) 

Cprv 

<!4> 

Chxe  ■ 

(46) 

For  three  dimensions  3D  acoustic  code  is  arguably  simpler  than  the  electromagnetic  case  since 
there  are  not  two  vector  fields.  However  porting  3D  electromagnetic  algorithms  to  the  acoustic 
case  is  not  as  trivial  as  in  two  dimensions. 

3  The  AFP  TFSF  Method 

Any  implementation  of  a  TFSF  boundary  requires  knowledge  of  the  incident  field  at  nodes  adjacent 
to  the  boundary  for  every  time-step  of  the  simulation.  Conceptually,  the  AFP  version  of  the  TFSF 
method  is  quite  simple  in  that  it  parallels  the  usual  description  of  propagation  in  the  continuous 
world. 

In  the  continuous  world,  the  spatial  dependence  of  a  harmonic  plane  wave  is  given  by  exp(— j\c- 
r)  where  k  is  the  wave  vector  and  r  is  the  position  vector  (exp (jut)  temporal  dependence  is 
understood).  For  a  pulsed  plane  wave  each  spectral  component  can  be  weighted  by  the  appropriate 
amount  to  give  the  frequency-domain  representation.  So,  for  example,  if  a  field  component  were 
found  at  the  origin  to  be  given  in  the  time-domain  by  /(f),  its  frequency-domain  representation 
would  be  F(u)  =  F[f(t)\  where  is  the  Fourier  transform.  The  field  at  an  arbitrary  point  r 
merely  has  to  account  for  the  displacement  from  the  origin.  Thus  in  the  frequency  domain  the  field 
is  given  by  F(r,  u)  =  F(u)  exp(—jk  •  r).  The  time-domain  signal  at  r  is  the  inverse  transform  of 
this,  i.e., 

f  (r,  t)  =  —  /  F(u)e~jkTejutdw.  (47) 

2^  J  —  OO 

Assuming  lossless  media,  in  the  continuous  world  the  magnitude  of  the  wave  vector  is  given  by  uj/c 
where  c  is  the  speed  of  light.  Thus  it  is  straightforward  to  evaluate  (47) — the  complex  exponential 
involving  space  merely  represents  a  shift  operation. 

In  the  discretized  world  of  the  FDTD  method,  one  can  follow  steps  which  parallel  those  in  the 
continuous  world.  However,  finding  the  field  at  an  arbitrary  point  is  complicated  by  the  fact  that 
the  wave  vector  in  the  grid  is  governed  by  the  FDTD  dispersion  relation  which  does  not,  in  general, 
have  a  closed-form  solution.  Nevertheless,  it  is  relatively  easy  to  solve  for  the  wave  vector  and  to 
perform  the  necessary  transforms  to  calculate  the  incident  field  wherever  it  is  needed. 

Since  the  input  signal  is  not  periodic,  one  cannot  use  a  discrete  Fourier  transform  (which  in¬ 
herently  assumes  a  periodic  signal).  Instead,  for  a  transient  signal,  one  must  use  the  discrete-time 
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Fourier  transform  [8].  The  inverse  and  forward  transforms  of  the  discrete  signal  f(qAt)  =  f[q]  are 
given  by 

M  =  jjjf  (48) 

oo 

F(<J)  =  £  f[q]e (49) 

g=—oo 

Note  that,  despite  this  transform  pertaining  to  a  discretized  world,  the  frequency  u>'  is  a  continuous 
variable.  In  FDTD,  any  “source  function”  f[q]  can  be  assumed  to  start  at  q  =  0  and  be  limited  to 
a  maximum  of  Ns  time  steps  (i.e.,  the  source  function  has  either  decayed  to  zero  at  time-step  Ns 
or  is  switched  off  at  time-step  Ns — this  could  correspond  to  the  time  at  which  the  simulation  is 
terminated).  Thus  (49)  can  be  written 

Ns- 1 

*V)  =  £  (50) 

g=0 


For  the  sake  of  illustration,  consider  a  ID  computational  domain  in  which  the  source-function 
f[q]  represents  the  time-varying  field  at  some  point.  Further  assume  this  field  is  propagating  in  the 
positive  x  direction.  The  field  at  a  point  which  is  m  spatial  steps  away  from  where  f[q]  is  measured 
would  be  given  by 

f[m,q}  =  ~  /  F{J)e-iHmejw'qdJ  (51) 

Jo 

where  k  is  the  numeric  wavenumber  and  5  is  the  spatial  step  size.  Equation  (51)  is  the  ID  FDTD 
analog  of  (47).  In  one  dimension  a  closed-form  expression  can  be  obtained  for  k5: 


(52) 


where  Sc  is  the  Courant  number.  (One  can  equate  J  with  the  more  familiar  uAt  which  typically 
appears  in  the  FDTD  dispersion  relation,  but  the  fact  remains  that  u /  varies  continuously  in  the 
integral  of  (51).) 

Knowing  the  direction  of  propagation  and  the  incident  field  at  a  given  point  (i.e.,  f[q]),  one 
can  calculate  F(u')  using  (50)  and  then  find  the  field  at  an  arbitrary  point  f[m,  q]  using  (51).  This 
was  done  in  [9]  where  it  was  shown  how  this  approach  could  predict  the  superluminal  component 
of  FDTD  propagation.  It  was  speculated  in  [9]  that  this  technique  could  be  used  to  realize  a 
TFSF  boundary  which  would  essentially  be  perfect.  Indeed,  this  approach  is  at  the  core  of  the 
work  by  [4]  and  [5],  but  neither  paper  provided  details  concerning  the  Fourier  transforms  nor  were 
details  discussed  in  [9].  (It  should  also  be  pointed  out  that  Ma  et  al.  [10]  have  obtained  an  analytic 
expression  for  the  field  at  an  arbitrary  point  in  the  FDTD  grid  due  to  a  source  which  is  impulsive 
both  in  time  and  space.  That  differs  from  the  solution  here  in  that  plane  waves  are  of  interest. 
These  plane  waves  may  be  impulsive  in  time  but  cannot,  by  definition,  be  impulsive  in  more  than 
one  direction.) 
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(53) 


Substituting  (49)  (after  a  change  of  index  from  q  to  q')  into  (51)  and  rearranging  yields 


f[m,q\ 


e-jk6meju>'(q-q')du,' 


This  equation  is  exact  and  the  integral  possesses  some  interesting  properties  (for  example,  it  is 
zero  when  m  >  (q  -  q')).  Note  that  the  source  function  f[q']  is  assumed  to  be  zero  for  q'  > 
Ns  but  f[m,  q]  can  be  evaluated  for  any  value  of  q,  i.e.,  it  is  not  bound  by  Ns.  Unfortunately 
this  equation  cannot  be  easily  evaluated  efficiently  nor  can  it  be  evaluated  without  resorting  to 
numerical  approximations. 

An  efficient  calculation  of  this  expression  is  obtained  by  employing  standard  discrete  Fourier 
techniques  (as  was  done  in  [4,5,9]).  This  is  equivalent  to  approximating  the  integral  in  (51) 
as  a  Riemann  sum.  Let  us  assume  that  Nt  equally  spaced  samples  of  the  integrand  are  used  to 
approximate  the  integral.  In  this  case  the  frequency  J  is  given  by  2irn/Nt  (where  0  <  n  <  Nt  - 1) 
and  du'  is  approximately  2n/Nt.  An  approximation  of  (51)  is  thus 


e-jk8mej^q 


2lX 

N't’ 


(54) 


Regrouping  terms  and  employing  (50)  yields 


(55) 


If  one  sets  Ns  equal  to  Nt  (which  can  be  accomplished  by  zero-padding  the  source  function  f[q']  to 
the  necessary  length),  then  the  term  within  braces  is  recognized  as  the  discrete  Fourier  transform 
of  the  source  function.  This  transform  is  multiplied  by  exp (—jk6m)  and  then  the  inverse  discrete 
Fourier  transform  is  taken.  Since  discrete  Fourier  transforms  can  be  calculated  efficiently,  (55)  can 
be  calculated  efficiently. 

In  the  implementation  of  a  TFSF  boundary,  one  must  calculate  fields  which  are  offset  both 
spatially  and  temporally.  Even  though  two  fields  are  temporally  offset  by  half  a  time  step  they 
still  use  identical  samples  of  the  source  function  f[q'].  The  offset  is  accounted  for  in  the  inverse 
transform,  i.e., 


r  n  Nt- i 

f  m  +  l,q  +  l  = 

n=0 


=  Y  F{un)e-j'kS/2ej^e-j'kSmej^ 


where  F(un)  is  the  term  in  braces  in  (55).  (When  calculating  the  various  field  components  one 
must  also  account  for  the  characteristic  impedance  and  the  orientation  of  the  fields.  This  is  dis¬ 
cussed  more  later.) 

One  question  which  remains  is  the  value  which  should  be  chosen  for  Nt  to  obtain  a  good 
approximation  of  the  exact  integral.  Naturally,  the  more  points  the  better  the  approximation  will 
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be,  but  one  can  obtain  rough  guidelines  as  follows.  Using  discrete  transforms  is  equivalent  to 
assuming  the  source  function  is  periodic.  However,  one  does  not  want  this  periodic  behavior  to  be 
evident  in  the  simulation.  Hence  the  discrete  Fourier  transforms  must  be  long  enough  (i.e.,  Nt  must 
be  great  enough)  so  that  virtually  all  the  energy  associated  with  the  incident  field  has  traversed  the 
total-field  region  before  the  incident  field  can  repeat  itself. 

For  example,  consider  an  incident  pulse  which  is  non-zero  for  100  time-steps  (Ns  =  100) 
and  which  is  incident  upon  a  TF  region  which  is  50  spatial  steps  wide.  The  inverse  transform 
associated  with  the  last  point  in  the  TF  region  must  be  sufficiently  long  so  that  it  can  model  the 
time  it  takes  for  the  incident  pulse  to  propagate  to  that  point  and  the  time  it  takes  for  the  pulse  to 
completely  pass  this  point.  Although  the  pulse  started  by  being  bound  by  100  time  steps,  because 
of  the  dispersion  in  the  FDTD  grid,  it  will  take  more  than  100  time  steps  for  the  pulse  to  pass  any 
point.  The  more  distance  the  pulse  has  to  travel,  the  more  it  will  disperse.  In  fact,  as  discussed 
in  [5, 11, 12],  the  group  velocity  in  the  FDTD  grid  goes  to  zero  at  the  coarsest  discretizations  and 
hence  it  arguably  takes  an  infinite  time  for  a  pulse  to  pass  completely  any  point.  Nevertheless,  in 
practical  applications  the  coarsest  discretization  are  not  of  interest.  By  using  a  reasonable  Courant 
number  and  a  reasonable  discretization  of  the  incident  pulse,  there  will  be  little  spectral  content 
at  the  coarsest  discretization.  Thus,  in  this  example,  a  discrete  Fourier  transforms  of  1024  points 
(i.e.,  Nt  =  1024)  would  almost  certainly  be  sufficient  to  describe  the  incident  pulse  over  the  entire 
TFSF  boundary. 

The  FDTD  simulation  itself  can  proceed  for  any  number  of  time  steps.  If  a  highly  resonant 
structure  were  being  illuminated  and  the  user  wanted  to  run  the  simulation  for  a  hundred-thousand 
time  steps,  or  more,  that  would  be  irrelevant  to  the  implementation  of  the  TFSF  boundary.  The 
incident  field  on  the  boundary  would  merely  be  assumed  to  be  zero  after  1024  time  steps. 

Generalizing  (51)  to  higher  dimensions  is  trivial  in  that  it  is  nearly  identical  to  (47) — the  only 
differences  are  the  limits  of  the  integral  and  the  use  of  the  discrete  wavenumber  components.  For 
example,  assuming  a  uniform  grid  in  which  Ax  =  Ay  =  <5,  let  the  field  f[m,  n,  q]  represent  the 
fields  at  the  point  r  =  (mS,  nS )  and  time  qAt.  Given  the  field  f[q]  at  the  origin  which  has  Fourier 
transform  F( u'),  f[m,  n,  q]  is  given  by 

1  r2n 

f[m,n,q]  =  —  /  F(u’)e~j*T  eju'qdu>' .  (57) 

27r  J  o 

In  this  case  the  components  of  the  numeric  wave  vector  k  must  be  calculated  from  the  the  2D 
dispersion  relation  but  the  integral  can  again  be  approximated  with  discrete  transforms  as  done  in 
(55).  We  note  that  the  superluminal  wave  vector  components  discussed  in  [9, 13]  are  not  incorpo¬ 
rated  in  the  results  to  be  shown  later.  These  components,  which  occur  at  the  coarsest  discretizations 
supported  by  the  grid,  experience  exponential  decay  as  they  propagate.  Discarding  them  from  the 
solution  slightly  increases  the  amount  of  leaked  fields  but  this  is  not  a  concern  in  practice  (owing  to 
the  associated  frequencies  not  being  ones  which  would  be  of  interest  and  the  inherent  exponential 
decay). 

Equation  (51)  gives  the  spatial  and  temporal  dependence  of  a  single  field  component.  Given 
a  single  field  component,  the  polarization,  and  the  direction  of  propagation,  all  the  other  field 
components  can  be  computed.  Despite  the  characteristic  impedance  of  the  FDTD  grid  being  exact, 
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such  a  computation  is  more  involved  than  in  the  continuous  world  because  of  the  non-orthogonality 
of  the  electric  field,  magnetic  field,  and  wave  vector.  Details  of  the  relationship  between  these  three 
quantities  are  discussed  in  [5]. 

When  a  halfspace  discontinuity  is  present,  one  must  account  for  the  reflected  or,  where  appli¬ 
cable,  the  transmitted  fields.  To  demonstrate  this,  we  first  consider  the  case  of  illumination  of  a 
Dirichlet  planar  boundary,  i.e.,  a  PEC  plane  where  there  is  no  transmitted  field  and  the  reflection 
coefficient  is  -1.  We  then  consider  penetrable  media  which  must  incorporate  the  transmission  and 
reflection  coefficients. 


4  TEZ  Polarization  and  a  PEC  Plane 

Consider  a  2D  FDTD  grid  with  TEZ  polarization  in  which  a  PEC  plane  is  assumed,  at  least  insofar 
as  the  incident  field  is  concerned,  to  span  the  computational  domain.  This  scenario  is  depicted 
in  Fig.  3.  The  acoustic  pressure  can  be  equated  with  the  magnetic  field  and  the  velocity  with  the 
electrics  field  (with  an  appropriate  change  of  signs).  The  PEC  plane  is  thus  equivalent  to  a  rigid 
boundary  at  which  the  normal  component  of  velocity  goes  to  zero.  We  define  the  “incident  field” 
as  being  the  sum  of  the  incoming  field  (i.e.,  the  field  whose  x -component  of  the  wave  vector  is 
positive)  and  the  reflected  field.  By  doing  this,  in  the  absence  of  any  additional  scatterers  other 
than  the  PEC  plane  itself,  the  SF  region  would  contain  no  fields.  It  is  important  to  note  that  in 
any  given  simulation  there  are  no  restrictions  on  the  contents  of  the  TF  region.  In  fact,  the  PEC 
plane  does  not  even  have  to  span  the  TF  region.  Thus,  for  example,  one  can  consider  the  fields 
associated  with  obliquely  illuminated  apertures.  In  the  SF  region,  however,  it  is  required  that  the 
plane  is  intact  and  that  no  other  scatterers  are  present. 

Nodes  which  are  tangential  to  the  TFSF  boundary  and  have  a  neighboring  node  on  the  other 
side  of  the  boundary  must  have  their  update-equations  corrected  to  account  for  the  existence  of  the 
TFSF  boundary  (see  [2]  for  details).  For  the  situation  considered  here,  the  TFSF  boundary  is  only 
three-sided.  The  field  is  not  specified  on,  or  beyond,  the  PEC.  Instead,  as  is  usual  when  modeling 
PEC’s,  the  tangential  electric  fields  along  the  PEC  are  set  to  zero  and  the  FDTD  algorithm  handles 
the  rest. 

To  implement  the  TFSF  method,  the  incident  field  contains  two  plane  waves:  the  “incoming 
field”  and  the  reflected  field.  This  scenario  is  shown  in  Fig.  4.  The  “reference  point”  in  Fig.  4  is 
essentially  the  origin  at  which  the  user  would  specify  the  source  function  f[q}.  Note  that  this  point 
does  not  need  to  be  within  the  computational  domain!  Its  location  is  on  the  PEC  and  such  that 
for  the  given  incident  angle  fa  at  the  start  of  the  simulation  the  incoming  field  has  not  yet  entered 
the  TF  region  (but  the  leading  edge  of  the  field  sits  poised  to  enter  the  region).  The  field  at  points 
along  the  TFSF  boundary  are  then  determined  relative  to  displacement  from  this  reference  point  in 
accordance  with  (57). 

One  can  easily  show  that  the  FDTD  reflection  coefficient  for  a  PEC  plane  is  identically  —1 
(relative  to  the  electric  field)  and  the  details  will  not  be  presented  here.  Suffice  it  to  say  that 
the  incoming  and  reflected  fields  have  the  same  y  components  of  their  wave  numbers  and  equal 
magnitude  but  opposite  signs  for  the  x  components.  The  Fourier  transform  of  the  source  function 
f[q\  is  taken  to  be  the  spectral  representation  of  the  magnetic  field  for  both  the  incoming  and 
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Figure  3:  Depiction  of  a  TEZ  grid  containing  a  PEC  plane.  The  TFSF  boundary  is  drawn  with  a 
dashed  lines.  The  nodes  enclosed  in  rounded  rectangles  must  have  their  values  corrected  owing 
to  a  neighboring  node  being  on  the  other  side  of  the  boundary.  In  the  acoustic  analog  of  this 
polarization  the  magnetic  field  can  be  equated  with  the  pressure  and  the  electric  field  with  the 
velocity.  The  PEC  plane  is  thus  equivalent  to  a  rigid  boundary. 
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reflected  fields.  The  x  and  y  components  of  the  electric  field  are  found  in  accordance  with  the 
orthogonality  condition  discussed  in  [5].  (The  incoming  and  reflected  fields  have  the  same  sign  for 
the  x  component  of  the  electric  field  but  opposite  signs  for  the  y  component.)  Thus  the  incident 
field  is  realized  by  summing  two  plane  waves,  each  with  a  common  reference  point  and  the  same 
amplitude.  The  only  difference  is  the  direction  of  propagation.  The  calculation  of  each  individual 
plane  wave  follows  the  details  provided  in  [5].  The  superposition  of  these  waves  satisfies  the 
boundary  condition  dictated  by  the  existence  of  a  PEC  plane  in  the  FDTD  grid  and  hence  can  be 
used  to  realize  a  nearly  perfect  TFSF  boundary. 

To  illustrate  the  behavior  of  this  TFSF  implementation,  consider  a  pulsed  plane  wave  propagat¬ 
ing  at  an  incident  angle  of  60  degrees.  The  field  is  traveling  in  free  space  and  the  Courant  number 
is  95  percent  of  the  3D  limit,  i.e.,  Sc  =  0.95/\/3.  This  Courant  number  was  chosen  so  that  the  re¬ 
sults  pertain  to  3D  simulations  in  which  the  incident  field  propagates  orthogonal  to  one  of  the  axes 
(the  factor  of  95  percent  was  chosen  somewhat  arbitrarily — stability  is  no  more  of  a  concern  with 
this  TFSF  technique  than  it  is  with  the  surrounding  grid).  The  pulse  is  a  Ricker  wavelet  discretized 
so  that  there  are  10  cells  per  wavelength  at  the  most  energetic  frequency.  With  this  discretization 
there  is  a  substantial  amount  of  energy  at  coarser  discretizations  (see  [5]  for  further  discussion 
of  this  source  function).  One  would  anticipate  this  pulse  would  suffer  substantial  dispersion  as  it 
propagates,  something  which  is  undesirable  in  practice  but  good  for  testing  the  performance  of  a 
TFSF  implementation.  Fig.  5  shows  a  computational  domain  which  is  180  by  200  cells.  A  vertical 
PEC  exists  105  cells  from  the  left.  The  SF  region  is  15  cells  thick.  Figure  5(a)  shows  the  magnetic 
field  150  time-steps  into  the  simulation.  The  images  show  the  log  base  10  of  the  absolute  value 
of  the  field  and  have  been  scaled  so  that  it  is  visible  over  three  orders  of  magnitude.  In  Fig.  5(a) 
the  incoming  field  has  already  encountered  the  PEC.  The  TFSF  boundary  is  aware  of  the  reflected 
field  and  there  is  virtually  no  leakage  through  the  boundary.  Fig.  5(b)  shows  the  magnetic  field  at 
350  time  steps.  The  dispersion  of  the  incoming  field  is  clearly  evident  as  the  width  of  the  pulse 
is  greater  than  it  had  been  at  150  time-steps.  This  dispersion  is  subsequently  evident  in  the  field 
reflected  by  the  PEC  plane.  The  AFP  TFSF  implementation  automatically  incorporates  these  nu¬ 
meric  artifacts.  Using  a  discrete  Fourier  transform  of  1024  points,  the  peak  magnetic  field  leaked 
across  the  boundary  in  this  case  is  approximately  five  orders  of  magnitude  down  from  the  peak 
value  of  the  magnetic  field  (i.e.,  100  dB  down  from  the  peak).  If  one  were  to  use  a  more  reasonable 
discretization  of  20  cells  per  wavelength  at  the  peak  frequency  of  the  pulse,  the  leaked  field  drops 
to  180  dB  down  from  the  peak  of  the  incident  field.  (Note  that  in  Figs.  5(a)  and  (b)  there  is  no 
reason  to  extend  the  computational  domain  beyond  the  PEC  boundary  since  no  fields  propagate 
past  the  PEC.  It  was  done  here  merely  for  the  sake  of  consistency  with  (c)  and  (d).) 

Figure  5(c)  also  shows  the  magnetic  field  over  the  computational  domain  at  time-step  350,  but 
in  this  case  there  are  two  slits  in  the  PEC  plane.  Each  slit  is  10  cells  wide  and  their  centers  are 
separated  by  40  cells.  The  field  scattered  back  to  the  left  of  the  PEC  as  well  as  the  field  which 
passes  through  the  slits  are  clearly  evident.  The  implementation  of  the  TFSF  boundary  is  oblivious 
to  the  actual  contents  of  the  TF  region  (or  to  anything  beyond  the  PEC  plane). 

The  diffraction  from  infinite  wedges  was  studied  using  the  FDTD  method  in  [14, 15].  In  that 
work  the  TFSF  boundary  passed  through  the  perfectly  matched  layer  (PML)  which  terminated  the 
grid.  For  a  field  originating  in  the  PML,  an  amplification  factor  had  to  be  found  to  compensate 
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Figure  4:  To  model  a  PEC  boundary,  the  “incident  field”  contains  both  an  incoming  and  a  reflected 
field.  The  width  of  the  TF  region  is  h#.  The  origin,  or  reference  point,  for  the  sake  of  calculating 
the  incident  field  is  a  distance  h#  cot  (fa)  from  the  bottom  of  the  TF  region.  This  ensures  that  the 
incoming  field,  which  is  specified  by  the  source  function  f[q],  is  completely  outside  of  the  TF 
region  at  the  start  of  the  simulation.  The  bending  of,  and  gap  in,  the  PEC  boundary  is  used  to 
emphasize  that  there  are  no  restrictions  on  the  contents  of  the  TF  region.  Inhomogeneities  can  be 
present  throughout  the  TF  region. 
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Figure  5:  Snapshots  of  the  magnetic  field  at  time-steps  (a)  150  and  (b)  350.  (c)  Magnetic  field 
at  time-step  350  when  two  slits  are  present  in  the  PEC  plane,  (d)  Magnetic  field  at  time-step  350 
showing  diffraction  from  the  edge  of  a  semi-infinite  plane.  The  incident  angle  is  60  degrees  with 
respect  to  horizontal. 


for  the  PML  loss.  Using  the  AFP  TFSF  technique,  it  would  not  be  necessary  to  have  an  incoming 
field  start  within  the  PML  owing  to  the  fact  that  the  AFP  TFSF  technique  already  includes  the 
reflected  field.  This  is  illustrated  in  Fig.  5(d)  which  shows  the  diffraction  from  a  “knife  edge.” 
This  snapshot  is  also  of  the  magnetic  field  and  at  time-step  350.  In  this  case  the  TFSF  boundary 
is  two-sided.  One-side,  which  is  drawn  vertically,  is  terminated  one  cell  before  the  top  edge  of  the 
computational  domain  where  a  second-order  Higdon  ABC  is  used.  The  grid  termination  would 
benefit  from  the  use  of  a  PML  as  described  in  [14, 15],  but  no  amplification  factors  would  be 
needed  and  the  reflected  field  would  be  completely  present  at  the  start  of  the  TFSF  boundary.  Thus 
the  reflected  field  would  not  have  to  build  up  from  the  start  of  the  PEC  plane  which  is  within  the 
PML  nor  would  it  suffer  the  corresponding  diffraction  at  that  leading  edge — there  is  no  “leading 
edge”  of  the  PEC  in  the  AFP  TFSF  implementation  for  the  incoming  field  to  encounter.  Instead  of 
the  knife-edge  shown  here,  wedges  could  be  studied  just  as  easily  (and,  as  will  be  more  clear  after 
considering  dielectric  boundaries,  the  technique  can  be  applied  to  penetrable  wedges  too). 

The  AFP  TFSF  technique  allows  the  calculation  of  the  incoming  and/or  the  reflected  field  at 
an  arbitrary  point.  The  technique  does  not  care  if  the  point  is  actually  on  the  TFSF  boundary  or 
even  if  it  is  within  the  grid.  Thus,  when  it  comes  to  recording,  for  example,  the  diffracted  field, 
the  observation  points  can  be  placed  anywhere  in  the  computational  domain.  One  can  subtract 
the  incident  or  reflected  field  from  the  recorded  field.  In  this  way,  one  does  not  have  to  restrict 
observation  points  to  the  scattered-field  region. 

There  are  other  problems  which  could  benefit  from  the  application  of  the  AFP  TFSF  boundary 
described  here.  For  example,  it  provides  an  alternative  way  to  study  the  scattering  from  randomly 
rough  surfaces  than  the  one  presented  in  [16].  In  [16]  the  incident  field  employed  a  Gaussian- 
tapered  plane  wave  as  described  in  [17],  The  Gaussian  taper  was  necessary  to  minimize  diffrac¬ 
tion  errors  which  would  be  present  if  an  obliquely  incident  plane  wave  were  to  encounter  a  finite 
surface — only  a  finite  amount  of  the  rough  surface  can  be  included  in  any  particular  simulation.  A 
Gaussian-tapered  plane  wave  is  not  a  true  solution  to  the  wave  equation.  Additionally,  the  taper  is 
such  that  a  very  large  computational  domain  must  be  used  to  ensure  the  fields  are  small  at  either 
end  of  the  tapered  wave.  On  the  other  hand,  using  the  TFSF  implementation  described  here,  the 
surface  is  effectively  infinite  (although  the  surface  roughness  must  be  contained  within  the  TF  re¬ 
gion).  The  surface  roughness  would  have  to  be  “turned  on”  (i.e.,  ramped  up  and  down  so  that  it 
met  the  edges  of  the  planar  surface  at  the  TFSF  boundary),  but  this  can  be  done  in  much  less  space 
than  the  Gaussian  tapering  of  the  incident  plane  wave.  Unlike  with  a  Gaussian-tapered  incident 
field,  in  the  AFP  TFSF  method  the  incident  field  is  a  solution  to  the  wave  equation.  Thus  the  AFP 
TFSF  boundary  has  the  potential  to  provide  much  more  efficient  and  more  accurate  solutions  to 
these  types  of  problems. 


5  TEZ  Polarization  with  a  Dielectric  Halfspace 

For  a  field  incident  on  a  penetrable  halfspace,  one  must  account  for  both  the  reflected  and  trans¬ 
mitted  fields  and  hence  must  know  the  reflection  and  transmission  coefficients.  Consider  a  plane 
wave  propagating  obliquely  in  the  2D  TEZ  grid  shown  in  Fig.  6.  The  non-zero  fields  are  Ex,  Ey, 
and  Hz  (this  corresponds  to  the  polarization  identified  as  TM  in  [4]).  Again,  the  acoustic  pressure 
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Figure  6:  Portion  of  the  TEZ  grid.  An  interface  exists  at  x  =  0.  To  the  left  of  the  interface 
the  permittivity  and  permeability  are  ei  and  /x1(  respectively,  and  to  the  right  they  are  e2  and  /x2- 
Along  the  interface  the  permittivity  is  ea.  Magnetic  fields  are  unambiguously  in  the  first  or  second 
medium  and  hence  always  use  either  Hi  or  /x2.  The  acoustic  analog  of  this  simulation  equates  the 
pressure  with  the  electric  field  and  the  velocity  with  the  magnetic  field. 


can  be  equated  with  the  magnetic  field  and  the  velocity  can  be  equated  with  the  magnetic  field.  For 
the  sake  of  simplicity,  assume  a  uniform  grid  where  Ax  =  Ay  =  8.  The  computational  domain 
consists  of  two  half-spaces  where  the  permittivity  and  permeability  are  ex  and  /xx,  respectively,  to 
the  left  of  the  interface  at  x  =  0  and  e2  and  /i2  to  the  right.  Throughout  the  following  a  subscript  1 
will  be  used  to  indicate  quantities  to  the  left  of  the  interface  and  a  subscript  2  will  indicate  quanti¬ 
ties  to  the  right.  The  Ey  nodes  along  the  interface  have  a  permittivity  of  ea,  the  value  of  which  is 
left  arbitrary  for  now. 

We  adopt  the  discrete  calculus  notation  described  in  [5]  (which  differs  slightly  from  that  use 
in  [4])  and  start  with  the  description  of  an  arbitrary  harmonic  wave.  The  magnetic  field  is  given  by 

H  =  azHz^azH0e-j{LmS+kyn6)ejwqAt 

=  azHoe-*Tej“qA\  (58) 

where,  for  propagation  at  an  angle  4>  relative  to  the  x  axis,  the  numeric  wave  vector  k  is 

k  =  f  kx ,  ky')  =  k  (cos  <f>,  sin  <f>) ,  (59) 
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uj  is  the  frequency,  q  is  the  temporal  index,  and  m  and  n  are  the  spatial  indices  in  the  x  and  y 

directions,  respectively.  The  spatial  dependence  is  given  by  exp(— jk  •  r)  where  r  =  (mS.  n8)  ( m 

and  n  are  not  restricted  to  integer  values  and  can  be  offset  by  appropriate  fractional  amounts  to 
account  for  the  staggering  of  the  grid).  The  corresponding  electric  field  is  given  by 

E  =  a XEX  +  ayEy  =  {axE0x  +  ayEoy)e'jt'reJaJ9A‘, 

=  Eoe~jfcV^A‘.  (60) 

The  vector  E0  and  scalar  H0  are  constants  for  a  given  frequency  (but  are  themselves  functions  of 
frequency).  A  tilde  indicates  a  numeric  quantity  while  a  caret  implies  a  quantity  is  in  the  frequency 
domain  and  may  be  complex. 

Let  the  shift  operator  shift  the  £-index  by  +1/2  where  £  E  {x,  y,  t}.  For  example, 

s+E  =  E0e_^x(m+V2)<+fcyn5)e.?ajgAt  =  e~^xS^2E.  (61) 


Conversely,  s ^  shifts  the  £-index  by  —1/2.  The  discrete  difference  operator  is  defined  as 


For  plane-wave  propagation,  the  discrete  difference  operators  can  be  represented  by  multiplicative 
functions.  When  £  is  either  x  or  y  one  obtains 


E  =  -jK$  E. 


(63) 


Similarly,  the  temporal  finite  difference  yields 

—  a  2  (  uj  A#  \  a  ^ 

<9tE  =  j-^  sin  ( ■ J  E  =  W E-  (64) 

The  difference  operators  acting  on  the  magnetic  field  yield  similar  results.  In  terms  of  these  oper¬ 
ators  the  dispersion  relationship  is  given  by 

«V  =  Kl  +  Kl  (65) 

Ignoring  the  shift  operators  which  are  common  to  both  sides,  for  the  two-dimensional  propa¬ 
gation  which  pertains  here,  the  FDTD  harmonic  form  of  Ampere’s  law  can  be  written 

jtteE  =  -jK  x  H  =  —j&xKyHz  +  j&yKxHz,  (66) 


where  K  =  (Kx,  Ky)  (see  [5]  for  further  details  including  the  shift  operators  which  have  been 
dropped).  Thus,  the  components  of  the  electric  field  are  related  to  the  magnetic  field  via 


Kyfj 

K*H 

n~eHz- 


(67) 

(68) 
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As  mentioned  previously,  knowing  one  field  component,  the  polarization,  and  the  direction  of 
propagation,  one  can  obtain  all  the  field  components.  Equations  (67)  and  (68)  demonstrate  this  is 
true. 

To  solve  for  the  reflection  and  transmission  coefficients,  one  must  obtain  two  independent 
equations  relating  them.  As  in  the  continuous  world,  the  phase  of  the  incoming,  reflected,  and 
transmitted  fields  must  match  along  the  interface.  This  dictates  that  the  angle  of  reflection  must 
equal  the  angle  of  incidence.  Assuming  a  unit  amplitude  incoming  wave,  the  incoming,  reflected, 


and  transmitted  magnetic  fields  can  be  written,  respectively, 

H4  =  aze~jiiiT, 

(69) 

IT  =  -aAe"^-r, 

(70) 

H4  =  a  zftee-*tr, 

(71) 

where 

ki 

=  axk\  cos  fa  +  ayki  sin  fa  =  (kix,  ky ) 

(72) 

kr 

=  ~axki  cos  fa  +  ayki  sin  fa  =  (-kix,  ky ), 

(73) 

k< 

=  a xk2  cos  fa  +  SLyk2  sin  fa  =  ( k2x ,  ky), 

(74) 

ki  and  k2  are  the  FDTD  wave  numbers  in  the  first  and  second  media,  respectively,  fa  is  the  incident 
angle,  fa  is  the  transmitted  angle,  and  f,e  and  Tte  are  the  reflection  and  transmission  coefficients, 
respectively.  The  temporal  dependence  which  is  common  to  all  terms  has  been  dropped.  Because 
of  the  phase  matching  which  must  exist  at  the  interface,  k\  sin  fa  must  equal  k2  sin  fa,  i.e.,  the  y 
component  of  the  wave  vector  is  the  same  throughout  the  grid. 

The  total  field  in  the  first  medium  is  the  sum  of  the  incoming  and  reflected  waves,  i.e., 

H,  =  a2(e-A'-f,.e-A-),  (75) 

E,  =  (76) 

The  total  field  in  the  second  medium  is  given  by  the  transmitted  field.  The  transmitted  electric  field 
is 

E4  =  — — -  x  a zTtee-^r.  (77) 

iZ€2 

The  vectors  K\  Kr,  and  K4  are  given  by  (Klx,Ky)  (~Klx,  Ky),  and  {K2x,  Ky),  respectively. 
Because  the  tangential  phase  is  the  same  for  all  the  fields,  it  is  also  true  that  Ky  is  the  same  for  all 
the  fields. 

Note  that  only  Ey  nodes  are  present  at  the  interface.  Nevertheless,  as  in  the  continuous  world, 
the  tangential  electric  field  must  match  at  the  interface,  i.e., 


Ei 


cc— 0 


=  ay-E4 


i=0 


(78) 
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Expanding  these  terms  yields 


TS  K  *  ~ 

sm(<j>j)y  ,  p  -j'fci  sin(&)j/ 

r>  e  '  r>  1  te& 

SZ61  iZ6i 


r$, 


Tte  e~-*k2  sin^dy_ 


Because  the  phase  must  match  along  the  boundary,  the  complex  exponential  can  be  eliminated. 
Canceling  Q  this  equation  can  be  written 

1  +  f,  =  (80) 

^2  X-lx 

Another  equation  relating  the  reflection  and  transmission  coefficients  can  be  obtained  from  the 
update-equation  for  the  electric-field  nodes  on  the  interface.  The  relevant  equation  is  the  y  com¬ 
ponent  of  Ampere’s  law  evaluated  at  x  =  0.  This  was  given  in  (66)  but  that  assumed  propagation 
in  a  homogeneous  space  which  is  no  longer  pertinent.  Instead,  we  explicitly  write  the  spatial  finite 
difference: 


J  6  q  E y 


=  ~  dxHz\ 

\X 

=  -\(  Hi  I 


\x=-S/2 


Unlike  before,  the  spatial  finite  difference  operator  dx  cannot  be  expressed  directly  in  terms  of  K’s 
since  the  difference  involves  fields  on  either  side  of  the  interface.  The  electric  field  in  (81)  can  be 
represented  as  either  the  transmitted  field  or  the  sum  of  the  incoming  and  reflected  fields — the  same 
result  will  ultimately  be  obtained.  Using  the  transmitted  field  for  the  electric  field  and  discarding 
common  phase  terms  yields 

jtea^T'e  =  (ftee-jK2*  -  [>'*  -  Ttee~jK'*] )  ,  (82) 

where 


k\  cos(4>i)S 

^lx  =  2  ) 

k2  cos  ((j)t)5 
*2*  =  - 2 - ' 

Combining  (80)  and  (82)  and  solving  for  the  reflection  coefficient  yields 


.  K^eiKlx  _  K^e-jK2x  _  j^-KlxK2x8 

p  _  62 _ 61 _ J  6162  _ 

Z&Le-jKi*  +  Ex z.e-jK2x  +  jJxl-K1xk2x8 


(S3) 

(84) 

(85) 


As  the  discretization  goes  to  zero,  the  third  term  in  the  numerator  and  denominator  goes  to  zero,  the 
complex  exponentials  approach  one,  and  the  Kx  s  approach  the  x  component  of  the  wavenumber 
in  their  respective  media.  Thus  this  expression  gives  the  continuous-world  reflection  coefficient  as 
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the  discretization  goes  to  zero.  It  is  interesting  to  note  that  this  is  true  even  though  at  this  point  we 
have  placed  no  restrictions  on  ea,  the  permittivity  used  for  the  nodes  along  the  interface  (although 
inherent  in  the  derivation  is  the  restriction  that  ea  cannot  be  a  pathological  value  such  as  zero). 
Using  (80)  in  (85),  the  transmission  coefficient  is  found  to  be 


_ 2^*  cos(/cla) _ 

+  Elx.e-jK2X  -f  j-e*-KlxK2x8' 

62  £l  J  Cl  £2 


(86) 


Using  Euler’s  formula  to  expand  the  complex  exponentials  and  writing  the  K’s  in  terms  of 
sines,  the  reflection  coefficient  can  be  written 


sinfc^) cos(kxi)  _  sinf^i) cos(rea;2) 

_ £2 _ £l _ 

sin(>Cg2)cos(Kxi)  .  sinf/Cji)  005(^3,2) 
£2  ei 


-jtt 


+j* 


where  the  imaginary  term  in  the  numerator  and  denominator  is 

sin(KIi)  sin(/cz2)  t  . 

N  =  — 1 - L  (2e0  -  d  -  62) 


(87) 


(88) 


Note  that  when  ea  is  the  average  of  the  permittivities  to  either  side  of  the  interface  R  is  zero.  In  this 
way  the  imaginary  part  of  the  reflection  and  transmission  coefficients  is  zero.  Since  the  continuous- 
world  reflection  and  transmission  coefficients  are  purely  real,  the  imaginary  part  is  an  inherent 
error.  The  expressions  above  effectively  constitute  a  proof  that  using  the  average  permittivity  for 
the  interface  nodes  is  optimum.  For  ea  =  (ei  +  e2)/2  the  reflection  and  transmission  coefficients 
become 


ei  sin (kx2)  cos(Ksi)  -  e2  sin(/cxl)  cqs(kz2) 
e'i  sin(Kl2)  cos(/czl)  +  e2  sin(«a:1)  cos(/cx2)  ’ 

_ 2e2  sin(ftgl)cos(/cai) _ 

ei  sin(Kx2)  cos(/cxx)  +  e2  sin(Kxl)  cos(/cx2) ' 


(89) 

(90) 


Despite  the  change  in  permeability,  these  equations  are  seemingly  independent  of  permeability. 
However,  the  permeabilities  dictate  the  wave  numbers  in  the  different  media  and  hence  the  perme¬ 
abilities  are  implicitly  contained  within  these  equations. 

Reflection  and  transmission  coefficients  for  this  polarization  were  provided  in  [4]  (ref.  (72)  and 
(73)  of  that  paper).  The  reflection  coefficient  in  that  work  involves  the  sum  of  28  terms  (these  terms 
involve  the  product  of  a  total  of  17  complex  exponentials  and  30  sine  functions)  and  was  obtained 
with  the  aid  of  a  computer  algebra  package.  The  complexity  of  that  expression  is,  it  must  be  noted, 
a  consequence  of  Moss  et  al.  considering  more  complex  media  than  that  assumed  here.  However, 
their  final  results  do  tend  to  obscure  the  simplicity  which  pertains  to  the  problem  of  interest  here. 
The  reflection  coefficients  which  pertain  to  PML  interfaces  have  also  been  studied  extensively. 
Derivations  of  the  numeric  PML  reflection  coefficient  can  be  found  in  [18-21]. 

With  the  reflection  and  transmission  coefficients  known,  the  incident  field  can  be  calculated  at 
an  arbitrary  point  given  the  source  function  f[q]  at  a  reference  point  and  the  direction  of  propaga¬ 
tion  of  the  incident  field.  The  location  of  the  reference  point  is  unchanged  from  that  depicted  in 


24 


Fig.  4.  Here  we  take  the  “incident  field”  to  mean  the  sum  of  the  incoming  and  reflected  field  if  a 
node  is  to  the  left  of  the  interface  and  the  transmitted  field  if  the  node  is  to  the  right.  Thus,  the 
time-domain  incident  magnetic  field  for  points  to  the  left  of  the  interface  are  given  by 


Hlz 


(91) 


where  F{u)  =  E[f[q)}.  As  discussed  in  connection  with  (56),  the  offsets  of  1/2  are  accounted  for 
in  the  inverse  transforms.  From  (67)  and  (68),  the  x  and  y  components  of  the  electric  field  are 


Eu 


Ely 


m  +  -,n,q  = 


-T 


-l 


(e-jli-r  _  f(ee-^"  r) 


m,n  +  -,q  = 


T 


-i 


Kx-^~  (e_A'r  +  rtee-^r) 


(92) 


(93) 


These  expressions,  as  well  as  the  corresponding  expression  for  the  transmitted  fields,  are  evaluated 
for  the  points  adjacent  to  the  TFSF  boundary. 

To  illustrate  the  behavior  of  the  TFSF  boundary,  Fig.  7  shows  two  snapshots  of  the  magnetic 
field  in  a  computational  domain  which  is  180  by  200  cells.  Free  space  is  to  the  left  and  extends 
over  90  cells.  To  the  right  is  a  dielectric  with  a  permittivity  of  4e0  (Ey  nodes  along  the  interface  use 
a  permittivity  of  2.5e0).  The  permeability  in  both  regions  is  that  of  free  space.  This  is  equivalent 
to  an  acoustic  simulation  in  which  the  sound  speed  in  the  second  medium  is  half  that  of  the  first. 
The  incoming  pulse  is  a  Ricker  wavelet  discretized  such  that  there  are  20  cells  per  free-space 
wavelength  at  the  most  energetic  frequency.  This  corresponds  to  10  cells  per  wavelength  in  the 
dielectric.  The  SF  region  is  15  cells  thick,  the  incident  angle  is  60  degrees,  and  the  Courant  number 
is  0.95/\/3. 

Figure  7(a)  shows  the  Hz  field  at  150  time  steps  when  the  leading  edge  of  the  pulse  has  first 
encountered  the  dielectric  at  the  center  of  the  bottom  of  the  figure.  Figure  7(b)  shows  the  field  at 
350  time  steps.  Now  the  reflected  and  transmitted  field  are  clearly  evident.  Since  no  scatterer  is 
present  in  this  simulation,  no  scattered  fields  are  visible  in  the  SF  region  (the  plot  uses  three  decades 
of  logarithmic  scaling).  One  can  clearly  see  the  refraction  of  the  transmitted  field.  Also,  owing  to 
the  higher  permittivity  of  the  second  medium,  the  transmitted  field  suffers  more  numeric  dispersion 
than  fields  in  the  first  medium.  (Dispersion  in  the  FDTD  grid  is  dictated  by  the  discretization  [2]. 
The  higher  permittivity  in  the  second  medium  results  in  shorter  wavelengths,  and  hence  coarser 
discretization  and  greater  dispersion,  than  in  the  first  medium.)  The  increased  dispersion  causes 
the  transmitted  pulse  to  broaden  noticeably  as  it  propagates — one  can  see  that  the  pulse  is  thinest 
at  the  interface.  One  may  ask  why,  at  any  given  time  step,  the  incoming  way  does  not  have  a 
similar  appearance  to  the  transmitted  field,  i.e.,  thinner  to  the  left  and  thicker  to  the  right?  The 
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TFSF  boundary 
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Figure  7:  Snapshots  of  the  magnetic  field  at  time-steps  (a)  150  and  (b)  350.  This  is  equivalent  to 
snapshots  of  the  pressure  in  a  fluid-fluid  problem  where  the  sound  speed  in  the  second  medium  is 
half  of  that  in  the  first  medium 

answer  is  that  the  incident  field  along  the  TFSF  boundary  is  exactly  matching  the  phase  speed  for 
all  the  spectral  components  for  the  particular  incident  angle,  i.e.,  all  spectral  components  have  the 
same  fa.  At  the  interface  between  free  space  and  the  dielectric,  boundary  conditions  dictate  the 
fields  must  be  continuous.  However,  the  free-space  phase  speeds  are  not  matched  to  the  phase 
speeds  in  the  dielectric  so  as  to  yield  a  single  transmitted  angle.  Since  the  phase  speeds  are  a 
function  of  the  frequency,  this  causes  the  depth-dependent  broadening  (or  thought  of  another  way, 
the  frequency-dependent  refraction  where  fa  is  a  function  of  frequency). 

As  with  the  PEC  simulation,  the  leaked  fields  are  approximately  100  dB  down  from  the  peak 
interior  fields.  When  a  more  reasonable  discretization  is  used  (i.e.,  there  is  not  significant  energy 
with  discretizations  less  than  10  cells  per  wavelength  in  the  second  medium),  the  peak  leaked  fields 
are  more  than  180  dB  down  from  the  peak  of  the  incident  field. 

6  TM*  Polarization  with  a  Dielectric  Halfspace 

Consider  the  TM2  grid  shown  in  Fig.  8.  For  this  polarization  the  electric  field  can  be  equated  with 
pressure  and  the  magnetic  field  can  be  equated  with  the  velocity.  The  interface  is  aligned  with 
Hx  and  Ez  nodes.  The  arrangement  and  indexing  of  nodes  is  consistent  with  the  TE2  grid  in  that 
both  grids  could  be  considered  slices  of  a  3D  grid  where  magnetic-field  nodes  are  centered  on  the 
faces  of  the  Yee  cube  while  electric-field  nodes  are  centered  on  the  edges.  For  this  polarization 
both  electric-  and  magnetic-field  nodes  lie  on  the  interface.  The  permittivity  and  permeability 
associated  with  these  nodes  is  ea  and  fxa,  respectively,  which  are  left  arbitrary  for  now.  We  again 
wish  to  find  the  reflection  and  transmission  coefficients.  The  incoming,  reflected,  and  transmitted 
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Figure  8:  Depiction  of  a  TMZ  grid  with  two  dielectric  halfspaces.  This  is  equivalent  to  an  acoustic 
simulation  where  the  electric  field  can  be  equated  with  pressure  and  the  magnetic  field  nodes  can 
be  equated  with  the  components  of  the  velocity  vector. 
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electric  field  are  given  by,  respectively, 

Ei  =  aze"^r,  (94) 

Er  =  azf(me-^r,  (95) 

E*  =  ajtute-***.  (96) 

where  f  ,m  and  Ttm  are  the  TMZ  reflection  and  transmission  coefficients,  respectively.  Temporal 
dependence  is  given  by  exp(juqAt)  and  is  common  to  all  terms  (hence  it  is  not  shown  explicitly). 
The  definitions  of  other  terms  are  as  before.  The  FDTD  form  of  Ampere’s  law  is 

jVLeEz  =  dxHy  —  dvHx.  (97) 

Because,  for  the  assumed  geometry,  the  phase  dependence  in  the  y  direction  is  the  same  throughout 
the  grid,  dy  can  be  replaced  with  -jKv.  For  any  of  the  fields,  the  FDTD  form  of  Faraday’s  relates 
the  electric  and  magnetic  fields  via 

-jftH  =  -jK  x  E  =  ~axjKyEz  +  a yjKxEz  (98) 

so  that  the  magnetic  field  components  are  given  by 

*•  =  <»> 

ft,  =  ~^EZ.  000) 


As  before,  Kx  has  equal  amplitude  but  opposite  sign  for  the  incoming  and  reflected  waves. 

Matching  the  2  component  of  the  electric  field  at  the  interface  dictates  that  the  sum  of  the 
incoming  and  reflected  fields  must  equal  the  transmitted  field  at  x  —  0.  Since  the  phase  must  be 
equal  along  the  boundary  this  reduces  to 


i  +  r,m  —  Ttm. 


(101) 


The  other  equation  relating  the  reflection  and  transmission  coefficients  is  obtained  from  Ampere’s 
law  applied  to  the  nodes  on  the  interface.  Using  (99)  in  (97)  and  rearranging  yields 

/  K ^  \ 

j  (  Qea  -  — ^  )  Ez  =  dxHy  (102) 


I  x=-5/2 


(102) 


(103) 


where  H\y  is  the  y  component  of  the  sum  of  the  incoming  and  reflected  waves  and  Hy  is  the  y 
component  of  the  transmitted  field.  After  using  (100)  to  express  Hy  in  terms  of  Ez  and  discarding 
common  phase  terms,  (103)  yields 


,s  (n  V.  -  K2, )  T„ 


-—K2, Tme~‘^  +  —  A'ZI  (>■•  -  f 
/i2  Ml  ' 


(104) 
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where  the  k’s  are  defined  in  (83)  and  (84). 

Using  (101)  and  (104)  to  solve  for  the  reflection  coefficient  yields 

f  ^K'*ejKlx  ~  f2K^~jK2x  ~  JS  ftW,  -  K*) 
%Klxe-iKl*  +  ^K2xe~^  +  jS  (02/raeo  -  If2) ' 


(105) 


This  can  be  compared  to  (45)  of  [4]  which  appears  to  have  a  typographical  error.  (The  exponent  of 
the  first  term  in  the  numerator  has  the  wrong  sign.  Also  the  term  Kz  [with  no  numeric  subscript] 
which  appears  in  that  expression  is  never  explicitly  defined.)  The  transmission  coefficient  is 


2^K1xcos(k1x) 

^Klxe-iK'*  +  ^K2xe~^  +  jS  (f 22/rae0  -  if2) ' 


(106) 


When,  as  was  used  in  the  TEZ  case,  ea  is  the  average  of  the  permittivity  to  either  side  of  the 
interface,  one  can  write 

<5(Q2/xaea  -  if2)  =  ^  (Q2/iaei  +  ftVae2  -  2 if2)  (107) 

Assuming  the  permeability  Ha  is  the  harmonic  mean  of  the  permeabilities  to  either  side  of  the 
interface,  i.e.,  fia  =  2/j1/i2/(/0  +  H2),  the  right-hand  side  of  (107)  can  be  written 


^  (QW 1  -  Kl)  +  -  K £))  . 


2  \fi  i+/r2v  '  y'  H 1+H2 

Using  the  dispersion  relation  (65),  this  becomes 


<5/x2 


■Kl  + 


8  Hi 


■K\x. 


Hi  +  H2  Hi  +  H2 
Employing  the  definition  of  the  if ’s  (63),  allows  us  to  write  this  imaginary  term  as 

4/^2  •  2/„  \  ,  4/^i  2 


<5(^i  +  H2) 


sin  («ix)  -f 


<5(^1  +  H2) 


sin  (k2i). 


(108) 
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(110) 


This  final  form  is  convenient  because,  assuming  ea  =  (ei  +  e2)/2and/ia  =  2/i1/r2/(/r1 +/x2),  when 
the  complex  exponential  in  (105)  and  (106)  are  expanded  it  is  clear  that  the  imaginary  part  of  that 
expansion  identically  cancels  the  term  shown  in  (1 10).  Thus,  the  resulting  expressions  are  purely 
real  and,  after  employing  the  double-angle  formula,  given  by 


H2  sin^Kis)  -  Hi  sin(2K2x) 
H2  sin(2ACia.)  -f  Hi  sin(2 k2x)  ’ 

_ 2 H2  sin(2^ix) _ 

Hi  sin(2/cix)  +  H2  sin(2 k2x)  ’ 


(111) 

(112) 


These  expressions  reduce  to  that  of  the  continuous  world  when  the  discretization  goes  to  zero. 
Since  the  exact  expressions  are  purely  real,  this  serves  as  proof  of  the  optimality  of  using  the 
arithmetic  mean  for  the  interface  permittivity  and  the  harmonic  mean  for  the  permeability. 

Implementation  of  an  AFP  TFSF  boundary  for  the  TMZ  polarization  also  yields  leaked  fields 
which  are  approximately  100  dB  down  from  the  peak  field  when  the  incoming  field  is  very  coarsely 
discretized.  Using  a  discretization  which  is  typical  of  actual  practice,  the  leaked  field  obtains  a 
maximum  of  approximately  -180  dB  relative  to  the  peak  of  the  incoming  field. 
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7  Lossy  Materials 

This  section  discusses  the  construction  of  a  TFSF  boundary  which  pertains  to  the  case  of  a  plane 
wave  (identified  as  the  incoming  wave)  obliquely  incident  on  a  lossy  halfspace.  Although  not 
strictly  required,  the  incoming  wave  is  assumed  to  travel  in  a  lossless  medium.  As  before,  the 
boundary  between  the  two  media  is  assumed  to  be  planar  insofar  as  the  calculation  of  the  incident 
field  is  concerned. 

As  has  been  discussed,  the  AFP  TFSF  boundary  requires  that  one  be  able  to  calculate  analyti¬ 
cally  the  incident  field  at  an  arbitrary  point.  This  necessitates  calculation  of  the  numeric  wavenum¬ 
bers  as  well  as  the  numeric  reflection  and  transmission  coefficients.  (Note  that  for  any  particular 
simulation  the  interface  between  the  media  need  not  be  planar  and  the  contents  of  the  TF  region 
are  arbitrary.  The  assumption  of  homogeneous  halfspaces  and  a  planar  interface  are  germane  only 
to  the  calculation  of  the  incident  field.)  We  start  by  describing  the  equations  which  govern  FDTD 
propagation  in  lossy  medium. 

Consider  a  harmonic  field  polarized  in  the  2  direction  propagating  in  the  xy  plane  of  an  FDTD 
grid,  i.e.,  TMZ  polarization.  The  electric  and  magnetic  fields  can  be  written 

a  ZEZ  =  gLzEoe-jfcxm5+kyn8)ejuqAt 

=  a.zE0e~^'re’u’qAt,  (113) 

H  =  axHx  +  ayHy  =  (axH0x  +  ayH^)e~^T 

=  H0e-^V"9A*.  (114) 


where  S  is  the  spatial  step  size  (assumed  uniform  for  the  sake  of  notational  simplicity),  At  is  the 
temporal  step  size,  (m,  n)  are  the  spatial  indices  in  the  x  and  y  directions,  respectively,  q  is  the 
temporal  index,  co  is  the  frequency,  k  =  a xkx  +  a yky  is  the  wave  vector,  r  =  a xmS  +  a yn6  is  the 
position  vector,  and  E0  is  an  arbitrary  amplitude.  A  tilde  is  used  to  indicate  a  numeric  quantity, 
i.e.,  one  whose  value  in  the  FDTD  grid  differs  from  that  in  the  continuous  world.  The  amplitudes 
H0x  and  H0y  are  dictated  by  the  impedance  of  the  grid. 

Discretizing  time,  Ampere’s  law  expanded  about  the  time-step  q  +  1/2  is 


V 


x  h^/2  =  cEg+1  ~  ,  gE«+1+E« 

At  2 


(115) 


where  the  superscript  indicates  the  time  step,  <7  is  the  conductivity,  e  is  the  permittivity,  and  time¬ 
averaging  is  used  to  obtain  the  conduction  current  at  the  necessary  temporal  location.  For  the  given 
harmonic  fields,  the  temporal  finite  difference  can  be  expressed  as 


whereas  the  time-average  term  can  be  written 


E?+i/2  —  AE9+1/2. 


(116) 


(117) 
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The  discretized  form  of  the  curl  operation  is  unchanged  from  that  presented  in  [5].  For  a  harmonic 
field,  the  finite  differences  in  the  x  and  y  directions,  identified  as  dx  and  dy,  respectively,  can  be 
represented  by  multiplicative  operators,  i.e., 


(118) 


where  £  6  (x,  y}.  Thus  Ampere’s  law  can  be  written 


j KxH()y  -f-  jKyH0x  —  (jcfi  T"  <xA)Eq  —  j 60  1  j  Eq. 

The  phase  term  eju}(-q+1^2',At  was  common  to  all  terms  and  hence  dropped. 

Defining  K  =  a XKX  +  a yKy,  the  discrete  form  of  Faraday’s  law  can  be  written 

—jK  x  E9  = 


(119) 


(120) 


where  y  is  the  permeability  (it  is  assumed  the  magnetic  conductivity  is  zero  although  this  is  not 
required).  This  yields  two  equations 


4°  =  ~^Eo, 

Hy  0  =  -~E0. 


(121) 


(122) 


Using  (121)  and  (122)  in  (119),  one  obtains  the  dispersion  equation  for  lossy  media 


Kl  +  Kl  =  n  v  [l  -  3^ 


(123) 


This  is  the  dispersion  relationship  for  lossy  material  which  was  previously  considered  in  [22]  and 
more  recently  studied  in  [23].  We  define  the  complex  permittivity  to  be 


'  .oA 

£  =  eL  ’ 

=  e  1  —  jLc  cot 


(124) 


(125) 


where  the  loss-factor  La  is  defined  as  aAt/ (2e),  the  real  quantity  e  is  given  by  e0er  and  e0  is  the 
permittivity  of  free  space.  Expressed  in  terms  of  the  underlying  functions  the  dispersion  relation  is 

sin 2(kx)  +  sin 2(kv)  =  ^  sin2  1  -  jL a  cot  ,  (126) 

where  is  k^5/2  and  £  E  { x,y }. 

In  the  construction  of  the  AFP  TFSF  boundary  the  frequency  and  the  wave  vector  components 
in  the  first  medium  can  be  easily  calculated.  Due  to  phase  matching  along  the  boundary,  the 
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tangential  component  of  the  wave  number  must  be  the  same  in  both  media,  i.e.,  ky2  =  kyi  when  the 
interface  corresponds  to  a  constant  x  plane.  We  ignore  superluminal  wave  numbers  and  hence  ky2 
is  purely  real.  (Throughout  the  analysis  we  ignore  superluminal  propagation  which  is  inherently 
present  but  of  little  practical  concern.  See  [9, 13]  for  further  details.)  Since  the  frequency  is  know, 
the  only  unknown  in  (126)  is  kx,  the  normal  component  of  the  wavenumber  which  is  complex  due 

A  A 

to  the  loss.  We  separate  the  real  and  imaginary  parts  of  kx  as  kx  —  k'x  +  jk "  or,  correspondingly, 

kx8/2  =  kx  =  k'x  +  jn".  Plugging  this  into  (126),  expanding  terms,  and  separating  real  and 
imaginary  parts  yields 


where 


cos(4)  cosh(/c")  =  C", 
sin«)  sinh«')  =  C", 


C' 


C" 


1  +  2 
Ht^-t 

~~S[ 


•If  \  Hr&r  .  1 

sm  (Ky)  -  —  sin 
Lc  sin  (ui  At) . 


(127) 

(128) 

(129) 

(130) 


Solving  (127)  and  (128)  for  k'x  and  k"  (and  using  physical  arguments  to  eliminate  non-physical 
solutions)  yields 


where 


-cos 


UV  ' 
v/8  C'. 


) 


U  =  (m  +  y/(M  -  2(7')  (M  +  2 C")) 1/2 , 

K  =  M  —  y/{M  -  2C')(M  +  2C'), 

M  =  1  +  C'2  +  C"2. 


(131) 

(132) 

(133) 

(134) 

(135) 


The  numeric  reflection  and  transmission  coefficients  for  TM*  polarization  were  presented  in 
(111)  and  (112)  for  lossless  isotropic  media  and  in  [4]  for  uniaxial  media.  We  assume  an  inter¬ 
face  which  is  aligned  with  the  Ez  and  Hx  nodes  as  shown  in  Fig.  9.  The  FDTD  reflection  and 
transmission  coefficients  which  were  derived  already  still  pertain  to  the  lossy  case — the  only  dif¬ 
ference  is  that  the  permittivity  and  wavenumber  in  the  second  medium  become  complex.  When 
the  electric-field  nodes  on  the  interface  use  the  arithmetic  average  of  the  values  to  either  side  for 
the  conductivity  and  the  real  part  of  the  permittivity  while  the  magnetic-field  nodes  on  the  inter¬ 
face  use  the  harmonic  mean  for  the  permeability  (i.e.,  /ia  =  2^lyu2/(/x1  +  ^2)),  the  reflection  and 
transmission  coefficients  are,  respectively, 


/x2  sin(2ftia;)  -  Hi  sin(2fc2x) 

fi2  sin(2«ix)  +  Hi  sm(2k2x)  ’ 

_ 2 /i2  sin(2/cix) _ 

Hi  sin(2ftlx)  +  H2  sin(2£2x) ' 


(136) 

(137) 
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Figure  9:  Depiction  of  a  TMZ  grid  with  two  dielectric  halfspaces.  The  interface,  corresponding 
to  x  =  0,  is  aligned  with  Ez  and  Hx  nodes.  The  second  medium  is  assumed  lossy.  Along  the 
interface  the  electric-field  nodes  use  the  arithmetic  average  of  the  permittivity  and  conductivity  to 
either  side  while  the  magnetic-field  nodes  use  the  harmonic  mean  of  the  permeabilities.  The  dashed 
line  represents  the  TFSF  boundary.  The  nodes  enclosed  in  rounded  rectangles  are  tangential  to  the 
TFSF  boundary  and  hence  must  have  their  updates  corrected  using  the  incident  field  associated 
with  the  other  node  in  the  box. 
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These  values  reduce  to  the  continuous-world  values  in  the  limit  as  the  discretization  goes  to  zero.  If 
one  does  not  use  the  average  material  properties  at  the  interface.  This  additional  term  can  cause  the 
agreement  between  the  FDTD  and  continuous-world  values  to  be  better  than  when  using  averaging. 
However,  this  improvement  only  exists  over  a  very  narrow  band  of  frequencies  and  the  agreement 
is  worse  at  all  other  frequencies.  Hence  we  continue  to  assume  that  the  mean  values  are  used  for 
the  interface  material  parameters. 

Although  the  dispersion  relationship  was  derived  in  terms  of  TMZ  polarization,  the  same  results 
pertain  to  TEZ  polarization,  i.e.,  (131)  and  (132)  still  pertain.  The  reflection  and  transmission  coef¬ 
ficients  derived  previously  are  still  pertinent  provided  one  allows  the  permittivity  and  wavenumber 
to  be  complex.  The  assumed  TEZ  grid  is  shown  in  Fig.  10.  If  the  arithmetic  mean  is  used  for 
the  conductivity  and  permittivity  on  the  boundary,  the  reflection  and  transmission  coefficients  are 
essentially  unchanged  from  before,  i.e., 


£i  sm(kx2)  costal)  -  e2  sin(/ca:i)  cos (ftgg) 
£i  sin(«x2)  cos(kxi)  +  e2  sin(«si)  cos  (kx2)  ’ 

_ 2e2  sin(Ka;i)cos(/czi) _ 

€1  sin(Kl2)  cos(kxi)  -I-  e2  sin(Kxl)  cos(kx2) ' 


(138) 

(139) 


Using  the  wavenumbers  as  well  as  the  reflection  and  transmission  coefficients  described  here,  the 
implementation  of  the  AFP  TFSF  boundary  then  follows  the  implementation  described  previously. 


8  Incidence  beyond  the  Critical  Angle 

By  allowing  the  normal  component  of  the  wave  number  in  the  second  medium  to  be  complex,  the 
capability  is  automatically  present  to  model  incoming  fields  which  are  incident  beyond  the  critical 
angle,  i.e.,  the  fields  are  evanescent  in  the  second  medium.  As  discussed  in  [13],  the  critical  angle 
in  the  FDTD  world  differs  from  that  in  the  continuous  world  and  is,  in  fact,  a  function  of  frequency. 
Nevertheless,  we  will  refer  to  the  critical  angle  as  if  it  were  a  constant.  (When  the  second  material 
is  lossy,  the  concept  of  a  critical  angle  is  nebulous  since  some  energy  is  lost  in  the  second  medium. 
The  code  developed  here  handles  both  lossless  and  lossy  material  but  this  section  will  be  limited 
to  lossless  material.) 

When  modeling  incidence  beyond  the  critical  angle,  one  must  keep  in  mind  the  somewhat  un¬ 
usual  behavior  of  the  incident  field.  The  geometry  assumed  here  is  of  an  infinite,  pulsed  incoming 
plane  wave  propagating  obliquely  toward  an  infinite  planar  interface  separating  two  half  spaces  in 
which  the  speed  of  propagation  in  the  second  medium  is  faster  than  in  the  first.  Roughly  speaking, 
fields  in  the  second  medium  can,  and  will,  move  tangentially  along  the  boundary  faster  than  they 
can  in  the  first  medium.  However,  to  satisfy  the  boundary  conditions  along  the  interface,  these 
faster  moving  fields  will  couple  energy  back  into  the  first  medium.  These  fields  will  be  in  advance 
of  the  incoming  wave.  These  “advanced  fields,”  despite  arriving  at  any  given  point  before  the 
incoming  wave,  are  causal.  A  good  discussion  of  these  fields  can  be  found  in  [24,25]. 

In  theory,  the  advanced  fields  are  non-zero  throughout  space  and  this  could  be  problematic 
for  implementation  of  a  TFSF  boundary  which  assumes  the  fields  are  initially  zero  throughout 
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Figure  10:  Depiction  of  a  TEZ  grid  with  two  dielectric  halfspaces,  the  second  of  which  is  lossy. 
The  interface,  corresponding  to  x  =  0,  is  aligned  with  Ey  nodes.  For  the  interface  nodes  the 
conductivity  and  permittivity  use  the  arithmetic  mean  of  the  values  to  either  side. 
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the  computational  domain.  However,  in  practice  the  advanced  field  are  small  and  can  be  made 
arbitrarily  small  by  delaying  the  incoming  wave.  Additionally,  if  there  is  loss  present  in  the  second 
medium,  this  serves  to  diminish  the  advanced  fields. 

Despite  the  existence  of  these  advanced  fields,  the  AFP  TFSF  implementation  is  oblivious 
to  them — the  same  code  can  be  used  for  all  incident  angles.  To  illustrate  the  case  of  incidence 
beyond  the  critical  angle,  Fig.  1 1  shows  the  magnitude  of  the  Hz  field  in  a  TEZ  simulation  where 
Hi  —  fi2  =  fM 3,  en  =  2,  er2  =  1,  the  incident  angle  is  60  degrees,  and  the  Courant  number 
Sc  is  l/\/2.  The  computational  domain  is  180  x  200  cells,  the  interface  runs  vertically  through 
the  middle  of  the  grid,  and  the  TFSF  boundary  is  offset  15  cells  from  the  edge  for  the  grid.  The 
incoming  wave  is  a  Ricker  wavelet  discretize  such  that  there  are  20  cells  per  wavelength  (in  the  first 
medium)  at  the  frequency  with  the  greatest  spectral  content.  The  Ricker  wavelet  is  broad-banded 
(see  [5]  for  further  details  of  this  source  function). 

Figure  1 1(a)  shows  the  field  at  time-step  180,  shortly  after  the  field  has  become  clearly  visible 
in  the  TF  region.  The  incoming  wave  has  unit  peak  amplitude  and  the  images  use  logarithmic 
scaling  so  that  fields  are  essentially  visible  over  three  decades  (i.e.,  fields  greater  than  10~3  are 
visible).  In  Fig.  1 1(a)  one  can  see  the  distinct  incoming  field  as  well  as  the  “haze”  associated  with 
the  field  which  arrives  in  advance  of  the  incoming  wave.  This  leading  field  exists  throughout  the 
computational  domain  but  falls  off  as  one  moves  away  from  the  incoming  wave. 

Figures  11(b)  and  (c)  show  the  field  at  time-steps  330  and  460,  respectively.  Note  that  the 
incident  field  is  not  visible  in  the  SF  region.  For  this  particular  simulation  the  peak  leaked  field 
is  less  than  6  x  10-5.  This  is  much  greater  than  the  leaked  field  found  in  the  typical  beneath- 
critical-angle  case  where  the  leaked  field  is  on  the  order  of  10~9  for  reasonable  discretizations. 
The  amount  of  leaked  field  can  be  reduced  further  by  delaying  the  incoming  wave  or  by  changing 
the  source  function  so  that  it  is  has  less  low-frequency  content.  (Low-frequency  energy,  with  its 
long  wavelengths,  falls  off  very  slowly.)  If  loss  can  be  added  to  the  second  medium,  this  can  also 
significantly  reduce  the  leaked  field. 

To  verify  that  the  observed  leaked  fields  were  a  consequence  of  the  inherent  nature  of  the 
incident  field  and  not  the  result  of  a  coding  error,  a  simulation  was  performed  in  which  the  AFP 
technique  was  used  to  calculate  the  initial  field  at  every  point  within  the  TF  region.  In  this  case  the 
leaked  field  dropped  to  approximately  10-15.  Thus  the  implementation  was  correct  and  the  leaked 
fields  are  a  consequence  of  the  AFP  calculation  yielding  a  particular  non-zero  value  over  the  entire 
TFSF  boundary  but  the  TF  region  is  initially  zero.  Unfortunately  this  type  of  initialization  is  not 
a  practical  way  to  lower  the  leaked  field  for  two  reasons.  First,  it  is  computationally  cumbersome 
and,  second,  it  presupposes  no  scatterers  are  present  in  the  TF  region.  (Simulations  which  lack  a 
scatterer  are  of  no  practical  interest.) 

Figure  11(d)  is  also  a  snapshot  at  time-step  460.  However,  to  illustrate  that  the  contents  of 
the  TF  region  can  be  arbitrary,  a  notch  has  been  placed  in  the  interface.  The  notch  is  20  x  20 
cells  where  the  second  medium  now  protrudes  into  the  first.  Because  of  this  notch  the  fields  in  the 
second  medium  are  no  longer  purely  evanescent.  One  can  see  the  field  scattered  by  the  notch  and 
how  it  has  passed  into  the  SF  region  at  this  particular  time  step. 
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Figure  12:  Depiction  of  3D  computational  domain  with  associated  2D  auxiliary  simulations. 

9  Three-Dimensional  Simulations 

At  present  the  AFP  TFSF  method  cannot  be  applied  to  three-dimensional  (3D)  halfspace  problems 
with  arbitrary  incidence  because  the  numeric  reflection  and  transmission  coefficients  have  not  been 
obtained  for  the  general  3D  case.  Were  these  available,  the  3D  formulation  of  the  AFP  TFSF 
technique  would  be  a  trivial  extension  to  the  two-dimensional  (2D)  case.  One  would  merely  have 
three  components  to  the  wave  vector  instead  of  two  and  evaluation  points  would  have  three  indices 
instead  of  two.  However,  precomputing  the  incident  field  over  the  TFSF  boundary  in  3D  simulation 
would  potentially  be  costly.  In  3D  the  TFSF  boundary  is  six-sided  (as  depicted  in  Fig.  12).  There 
are  four  field  components  per  face.  If  the  incident  field  is  non-zero  over  a  large  number  of  time 
steps,  the  memory  required  to  store  the  incident  field  could  easily  exceed  the  memory  required  to 
store  the  fields  within  the  grid  itself.  Therefore  it  is  envisioned  that  a  general  3D  implementation 
using  the  AFP  technique  would  store  the  precomputed  field  to  disk  (possibly  stored  in  six  separate 
files  representing  the  six  faces).  In  this  case,  in  order  to  make  the  necessary  corrections  to  the  field 
tangential  to  the  TFSF  boundary,  a  3D  FDTD  simulation  would  merely  read  the  stored  incident 
field  from  disk.  This  shifts  the  calculation  burden  to  an  a  priori  step  and  the  memory  burden  to 
disk.  (Versions  of  the  authors’  2D  programs  allow  the  incident  field  to  be  precomputed,  stored  to 
disk,  and  then  read  concurrently  with  the  FDTD  simulation.) 

The  implementation  of  the  general  3D  problem  is  the  subject  of  on-going  research.  However, 
for  3D  problems  in  which  the  incident  field  propagates  orthogonal  to  one  of  the  grid  axes,  there 
is  a  relatively  simple  way  to  implement  the  AFP  TFSF  technique  efficiently.  The  approach  is 
illustrated  in  Fig.  12.  Here  it  is  assumed  that  the  incident  field  is  propagating  obliquely  to  the  x  and 
y  directions,  but  orthogonal  to  the  z  direction  (one  can  easily  permute  the  indices  for  propagation 
in  other  direction).  The  polarization  of  the  incident  field  is  arbitrary,  i.e.,  all  six  field  components 
are  allowed  to  be  nonzero.  Two  2D  auxiliary  simulations  are  performed  to  calculate  the  incident 
field,  one  for  the  TEZ  components  of  the  incident  field  and  one  for  the  TAP  components.  These 
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simulations  require  that  the  incident  field  be  specified  over  eight  one-dimensional  boundaries  (four 
for  each  of  the  2D  simulations). 

It  is  required  that  the  TF  regions  of  the  2D  simulations  be  at  least  as  big  as  the  cross-sectional 
area  of  the  TF  region  in  the  3D  grid  (with  the  geometry  shown  here,  this  would  correspond  to  the 
area  of  the  TF  region  seen  along  a  constant  z  plane).  The  auxiliary  simulations  do  not  include  any 
scatterers.  Since  the  fields  leaked  into  the  SF  region  by  the  AFP  technique  are  typically  so  small, 
one  need  not  pay  much  attention  to  termination  of  these  auxiliary  grids  (i.e.,  since  the  leaked  fields 
are  often  on  the  order  of  10-9  one  does  not  have  to  worry  about  using  ABC’s  to  terminate  the 
grids).  Note  that  the  assumed  locations  of  the  interfaces  shown  in  Fig.  9  and  10  are  consistent  with 
slices  through  a  3D  grid. 

Referring  to  the  geometry  shown  in  Fig.  12,  the  interior  of  the  TF  regions  of  the  auxiliary  grids 
would  provide  the  x  and  y  components  of  the  fields  needed  at  the  “top”  and  “bottom”  of  the  TFSF 
boundary  in  the  3D  grid  (i.e.,  at  the  two  constant  z  planes).  The  fields  on  the  constant  x  and  y 
planes  are  also  taken  from  the  auxiliary  simulations,  but  in  this  case,  since  there  is  no  variation  in 
the  z  direction,  the  value  is  only  a  function  of  x  or  y  (the  same  value  pertains  to  the  entire  vertical 
extent  of  the  face  of  the  boundary). 

In  the  FDTD  grid  the  decomposition  of  the  incident  field  into  TE  and  TM  components  is  not 
trivial  since,  in  general,  the  electric  field,  magnetic  field,  and  direction  of  propagation  do  not 
form  an  orthogonal  set.  One  can  specify  the  direction  of  travel  and  the  orientation  of  the  electric 
field,  but  then  one  cannot  use  a  simple  vector  projection  of  the  field  components  and  characteristic 
impedance  of  the  continuous-world  medium  to  find  the  relative  amplitudes  of  the  TE  and  TM  fields. 
Instead,  one  must  use  the  FDTD  impedance  relationship  to  determine  the  amplitudes.  Details 
concerning  polarization  can  be  found  in  [5]. 

10  Conclusion 

The  AFP  TFSF  method  can  be  readily  applied  to  problems  involving  planar  interfaces,  whether 
dielectric  or  PEC.  The  method  is  ideally  suited  to  oblique  incidence  and  does  not  suffer  the  in¬ 
herent  approximations  associated  with  using  an  auxiliary  grid.  If  no  wavenumber  components  are 
discarded  from  the  transforms,  the  only  errors  associated  with  the  technique  would  be  those  associ¬ 
ated  with  implementation  of  the  discrete-time  Fourier  transform.  In  practice,  even  when  discarding 
superluminal  wavenumber  components  and  using  a  coarsely  discretize  incoming  field,  the  leaked 
fields  are  approximately  100  dB  down  from  the  peak  excitation.  Using  a  discretization  which  is 
more  typical  of  actual  practice,  the  leaked  fields  are  approximately  180  dB  down.  Modularized 
programs,  written  in  C,  which  implement  the  AFP  TFSF  method  for  all  the  cases  considered  here 
are  available  from  the  PI. 

The  approach  used  here  is  not  restricted  to  the  second-order  Yee  FDTD  algorithm.  The  same 
steps  could  be  followed  to  derive  an  AFP  TFSF  technique  for  any  FDTD  method  which  has  a  rig¬ 
orous  dispersion  relation.  Additionally,  the  method  could  be  applied  to  multiple  layers  where  one 
would  have  to  solve  for  the  fields  in  the  multi-layer  system  as  is  done  in  the  continuous  world.  As 
will  be  discussed  in  a  companion  paper,  the  method  can  be  used  with  lossy  media,  can  model  inci¬ 
dent  angles  beyond  the  critical  angle  (i.e.,  where  the  fields  in  the  second  medium  are  evanescent), 
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and  can  be  implemented  efficiently  in  three  dimensions. 

For  problems  which  can  be  solved  both  by  the  AFP  technique  and  the  traditional  one-dimensional 
.  auxiliary-grid  approach,  the  AFP  technique  yields  far  greater  accuracy  (except  in  the  case  of  grid- 
aligned  propagation).  Furthermore,  the  AFP  TFSF  technique  provides  the  ability  to  study  many 
problems  which  cannot  be  solved  using  the  traditional  one-dimensional  auxiliary-grid  approach. 
For  incident  angles  beyond  the  critical  angle,  the  solution  obtained  from  the  AFP  technique  is 
compromised  somewhat  by  the  inherent  nature  of  the  field  which  arrives  in  advance  of  the  in¬ 
coming  wave.  This  degradation  is  unavoidable  given  the  seemingly  acausal  incident  field  and  it 
is  believed  that  no  other  TFSF  method  could  provide  better  fidelity.  The  AFP  technique  can  be 
applied  efficiently  to  3D  problems  in  which  propagation  is  orthogonal  to  one  of  the  axes.  For  elec¬ 
tromagnetics,  this  requires  that  two  auxiliary  2D  simulations  be  performed — only  one  would  be 
required  in  an  acoustic  simulation.  A  general  3D  implementation  appears  possible  but  this  awaits 
the  derivation  of  the  FDTD  reflection  and  transmission  coefficients  in  3D. 
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